2nd order initial value problem in matlab

[epic325]

Homework Statement

1.)I want to write a function in MATLAB that contains the 2nd order function:

20*d^{2}x;(dt^{2})+5*dx/dt + 20*x=0 (dampened spring)

-The function should have 2 inputs (time,[initial values]) initial values should be a vector of 2 values
-The function should output a vector of the time derivatives

2.)I want to write a script that that uses the function to plot the function using euler's method with a for loop.

Homework Equations

I solved for the time derivatives

[w(2) ; -.25w(2)-w(1)] where my initial values are [1; 0] x(0)=1 and dx/dt(0)=0

The Attempt at a Solution

I think I have figured the code out but it's not plotting the right thing. it should be a dampened oscillator, but its plotting a line. I just need some help with the debugging. I will paste the function code followed by the separate script:

%This is the function where the input will be t and wVec=[initial values]

function wDerVec = odeFunc(t,wVec)
wDerVec=[wVec(2) -.25*wVec(2)-wVec(1)];
end

clear
clc
h=0.1; %h is the time step (.1s)
t=0:h:15; %time interval (0 to 15 s)
wVec=[1 0]; %initial condition

for i=2:length(t),
odeFunc(t,wVec);
wVec(i,:)=wVec(i-1,:)+ans*h;
end

plot(t,wVec(:,1)); %plots the distance with respect to time
title('Euler Approximation, damped spring');
xlabel('Time');
ylabel('distance, x');


I'm not very good at MATLAB so any help would be greatly appreciated! Thanks!

 

[jfgobin]

Good Afternoon epic325,

I tried your equation in Octave and had no problem. Here is the sequence of commands I used.

function Y=odeFCT(x,t)
Y=zeros(2,1);
Y(1)=x(2);
Y(2)=(-1/20)*(5*x(2)+20*x(1));
endfunction
t=linspace(0,50,1001);
y=lsode(@odeFCT,[1 0]',t);
plot(t,y(:,1))

Here is the graph obtained (plotted with gnuplot instead of Octave.)

Let me dive into your code.

 

[jfgobin]

Quick question -

In your code, you have

for i=2:length(t),
odeFunc(t,wVec);
wVec(i,:)=wVec(i-1,:)+ans*h;
end

Shouldn't it be

for i=2:length(t),
odeFunc(t,wVec(i-1,:));
wVec(i,:)=wVec(i-1,:)+ans*h;
end

?

 

[epic325]

Yes, jfgobin. That was the problem. Thank you very much for your help!