Troubleshooting Simple Harmonic Oscillator in MATLAB

[ShayanJ]

I'm trying to plot the evolution of a simple harmonic oscillator using MATLAB but I'm getting non-sense result and I have no idea what's wrong!
Here's my code:

clear
clc
x(1)=0;
v(1)=10;
h=.001;
k=100;
m=.1;
t=[0:h:10];
n=length(t);

for i=2:n

	F(i-1)=-k*x(i-1);
	v(i)=(h/m)*F(i-1)+v(i-1);
	x(i)=h*v(i-1)+x(i-1);

end

plot(t(1:n),x);

And this is what I get:
http://www.cvberry.com/octave/graphs/graph1395294146456.jpg

What's wrong?
Thanks

 

[maajdl]

That's not very surprising.
You have discretized the differential equation.
You might solve the discretized version you got, and see that indeed its solution is what you found.

This means that you need to find a better numerical approximation.
One that is stable, or more stable.
This means an approximation that will be closer to the solution of the DE on a longer time span.

There is of course no approximation that will fit exactly the DE for any time.

There are many methods that improve your integration scheme.
One of the best known is the Runge-Kutta (RG) method.

In my personal toolbox methods are those two that I like to play (RG or other for serious work) with:

- "more precise stepping" (small improvement, less divergent)
v(i)=(h/m)*F(i-1)+v(i-1);
x(i)=F(i-1)/m*h²/2+v(i-1)*h+x(i-1);

- "implicit scheme" (decreasing amplitude)
by solving at each time step the linear system:
F(i)=-k*x(i);
v(i)=(h/m)*F(i)+v(i-1);
x(i)=h*v(i)+x(i-1);
(many variants are possible)

or (if I am not mistaken):

F(i) = -((m*(dt*k*v(i-1) + k*x(i-1)))/(dt**2*k + m))
v(i) = -((-(m*v(i-1)) + dt*k*x(i-1))/(dt**2*k + m))
x(i) = (m*(dt*v(i-1) + x(i-1)))/(dt**2*k + m)

- try RG: http://en.wikipedia.org/wiki/Runge–Kutta_methods

The stability of any numerical scheme can be analyzed.
You can also do numerical experiments, very easily with matlab.
You will see the benefit of renowned methods, like the RG.

- see also:
http://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations

- the DE can also be written as a 2x2 system of 1st order differential equations:

x' = v
v' = F/m = -k x / m

which can be put into matricial for:

transpose({x',v'}) = {{0,1},{-k/m,0}} transpose({x,v})

or: X' = M X where X = transpose({x',v'})

the solution of which is simply X = Exp(M t) Xo

Calculating Exp(M t) is almsot a standard of linear algebra,
and in the case where M is constant, you get an exact solution.
This is of course almost trivial here, as this would be equivalent to an analytical solution.
However, if M is not constant in time, then this integration scheme on a step by step basis is very efficient.
Also, Exp(M t) can also be calculated approximately, and a 1rst order Taylor gives you the method you already used.
Better approximations can for example be based on Padé approximants.At this stage, just try a Runge-Kutta method.

 

[ShayanJ]

I used RK and it worked,Thanks man!
Here's the new code:

clear
clc
x(1)=0;
v(1)=10;
h=.001;
k=10;
m=.1;
t=[0:h:10];
n=length(t);

for i=2:n-1

	k1=-((k*h)/m)*x(i-1);
	l1=h*v(i-1);
	k2=-((k*h)/m)*(x(i-1)+l1/2);
	l2=h*(v(i-1)+k1/2);
	k3=-((k*h)/m)*(x(i-1)+l2/2);
	l3=h*(v(i-1)+k3/2);
	k4=-((k*h)/m)*(x(i-1)+l3);
	l4=h*(v(i-1)+k3);
	
	x(i)=x(i-1)+(1/6)*(l1+2*l2+2*l3+l4);
	v(i)=v(i-1)+(1/6)*(k1+2*k2+2*k3+k4);

end

plot(t(1:n-1),x);

And the result:

http://www.cvberry.com/octave/graphs/graph1395303464048.jpg

I should confess I didn't think these methods are that much different!

 

[maajdl]

Click once on the "Thanks" button, just to increase my poor count a little bit.
Nice that you tried the RK!
But why does it work so better?
And on how many periods would it work so well?
150, 1500, 15000 ?
And why does the "simple" approach work so badly?
Any idea?

 

[ShayanJ]

Click once on the "Thanks" button, just to increase my poor count a little bit.
Nice that you tried the RK!
But why does it work so better?
And on how many periods would it work so well?
150, 1500, 15000 ?
And why does the "simple" approach work so badly?
Any idea?

I don't know why is that. I'm surprised too.
I plotted it till 500,it was OK.But this new code takes so much longer time to execute that I don't want to try larger values.
I'm just introduced to the numerical computation's world so I don't know enough to talk about it.

 

[maajdl]

To be honest, I don't know either.

 

[DrClaude]

The scheme you first use is called Euler's method. It is well known to be most often unstable: the inevitable numerical error that is made at each step grows exponentially, except for a small region of parameters.

 

[DrClaude]

But this new code takes so much longer time to execute that I don't want to try larger values.

There is a built-in ODE solver in MATLAB that will be much faster. Look up for instance the help page for rk45.

Also, MATLAB is notoriously slow when loops are used, if you are not careful. In your code, for instance, the vectors x and v change size every step of the loop (they increase in length by one). It is best to first initialize them to the right length (n-1, in your case) by using, for instance,

x = zeros(n-1,1);
v = zeros(n-1,1);
v(1) = 10;