CODEWITHSUNDEEP
matlabodenumerical-integration
whatsherface
whatsherface

Reputation: 413

Matlab: ode45 output incorrect for forced spring mass damper

Here is my Matlab code to solve the second order ODE for a mass-spring-dashpot system:

function Spring

clear all;
close all;
options=odeset('RelTol',1e-6); 
p0 = [1 0]; %initial position and velocity
[t,p] = ode45(@SpringFunction, [0 20], p0, options);

plot(t,p(:,1)) %plot the first column (x) vs. time

return

function pdot = SpringFunction(t,p)

c = 5; w = 2;
g = sin(t); % forcing function
pdot = zeros(size(p));
pdot(1) = p(2);
pdot(2) = g-c*p(2)-(w^2)*p(1);

return

I believe the output I am getting is wrong because for this case I think a plot of displacement vs. time should look like a sinusoid with decreasing amplitude. Instead, it looks like this: decreasing function which relaxes to a sinusoid whose amplitude is constant

This seems incorrect to me, but please correct me if I'm wrong. I cannot see what is incorrect in the code.

Upvotes: 3

Views: 2078

Answers (1)

horchler
horchler

Reputation: 18484

You're sinusoidally forcing a damped system so you should expect the steady state to be a sinusoid. Here's a nice vibrations tutorial (PDF) – see pp. 448–450 on damped sinusoidal forcing. Try removing the forcing or greatly reducing its amplitude. Also, it looks like you have a lot of damping. Your damping ratio, ζ (zeta), appears to be c/(2*w) = 5/4. This implies that the unforced form of your system is over-damped – with no forcing you won't see oscillation.

Also, be careful using ode45 with oscillatory systems. If your system happens to be too stiff, you may either need to adjust the tolerances or use a stiff solver such as ode15s to get accurate results. You can always try using tighter tolerances to check that they yield qualitatively similar results (same period, amplitude, rate of growth/decay, etc.).

Upvotes: 3

Related Questions