2nd Order ODE Euler method not producing expected results

Hi,

I am trying to solve a 2 degrees of freedom damped mass spring system using Euler's Method. I believe i am about 90% of the way there but the plot i am getting out is showing the displacement oscillations increasing exponentially(?) over time. I am pretty sure that the opposite should be happening.

Here is the starting matrix from which i've made my equations:

V and Z are respectively in order to solve with basic euler.

Here is my code:

close all;
clear;
clc
%Coefficients
m1=3;
m2=6;

k1=1;
k2=9;
k3=3;

c1=0.02;
c2=0.09;
c3=0.06;

% Step sizes
h=0.1;
t=0:h:100;
n=size(t);
N=n(2);

% Allocate array sizes by n steps
x1=zeros(n);
x2=zeros(n);
V=zeros(n);
Z=zeros(n);

% Initial Conditions
x1(1)=-7;
x2(1)=7;
V(1)=-3;
Z(1)=3;

% Create differential equations
dVdt=@(x1,x2,V,Z) 1/m1*(-(k1+k2)*x1+k2*x2-(c1+c2)*V+c2*Z);
dZdt=@(x1,x2,V,Z) 1/m2*(k2*x1-(k2+k3)*x2+c2*V-(c2+c3)*Z);



for i=1:N-1
    V(i+1)=V(i)+h*dVdt(x1(i),x2(i),V(i),Z(i));
    x1(i+1)=x1(i)+h*V(i);
    Z(i+1)=Z(i)+h*dZdt(x1(i),x2(i),V(i),Z(i));
    x2(i+1)=x2(i)+h*Z(i);

end

plot(t,x1, 'blue',t,x2,'red');
xlabel('t');
ylabel('x');

With the resulting plot ^^^

As you can see the x values become very large which definitely shouldn't be the case.

Can someone please point out the dumb thing I'm doing that's giving me this kind of result?

Thanks in advance!

Answer 1

[John D'Errico]

That is not uncommon with Euler. It can be an unstable method. At each step, a small amount of error gets added in. Eventually, the noise starts to dominate, if the step size is sufficiently large.

Better methods, (for example backwards Euler) were designed to avoid this problem.

One clue to this is to reduce the step size (h), and see if suddenly, the Euler method solution becomes well-behaved. For example, if I use h=0.0001 or 0.001, the result below is consistent. But you had h=0.1. This is apparently sufficiently large to push Euler into its unstable domain for this problem. Do some reading here:

https://en.wikipedia.org/wiki/Euler_method

m1=3;
m2=6;
k1=1;
k2=9;
k3=3;
c1=0.02;
c2=0.09;
c3=0.06;
% Step sizes
h=0.0001;
t=0:h:100;
n=size(t);
N=n(2);
% Allocate array sizes by n steps
x1=zeros(n);
x2=zeros(n);
V=zeros(n);
Z=zeros(n);
% Initial Conditions
x1(1)=-7;
x2(1)=7;
V(1)=-3;
Z(1)=3;
% Create differential equations
dVdt=@(x1,x2,V,Z) 1/m1*(-(k1+k2)*x1+k2*x2-(c1+c2)*V+c2*Z);
dZdt=@(x1,x2,V,Z) 1/m2*(k2*x1-(k2+k3)*x2+c2*V-(c2+c3)*Z);
for i=1:N-1
    V(i+1)=V(i)+h*dVdt(x1(i),x2(i),V(i),Z(i));
    x1(i+1)=x1(i)+h*V(i);
    Z(i+1)=Z(i)+h*dZdt(x1(i),x2(i),V(i),Z(i));
    x2(i+1)=x2(i)+h*Z(i);
end
plot(t,x1, 'blue',t,x2,'red');
xlabel('t');
ylabel('x');

My guess is this is probably what you were expected to trip over, why this problem was assigned to you. In general, Euler's method can be dangerous, certainly if you don't know what to watch out for.

[James Webb]

Thanks very much, this is what I suspected since the problem explicitly says to use a step size of 0.1 for euler, rk4 and ab4.

[John D'Errico]

Yep. They wanted you to have exactly this problem. I assume this would be a prelude to explaining exactly what we have told you, that forward Euler can be unstable when the step size is too large. The comparison to other solvers would show something strange, then making you endlessly search for what you did wrong. Thank your teacher. ;-)

This link does give a decent explanation of what happens, and why Euler can be unstable for some step sizes:

http://www.math.iit.edu/~fass/478578_Chapter_4.pdf

Answer 2

[Torsten]

Not that bad.

%Coefficients
m1=3;
m2=6;

k1=1;
k2=9;
k3=3;

c1=0.02;
c2=0.09;
c3=0.06;

% Step sizes
h=0.00005;
t=0:h:100;
n=size(t);
N=n(2);

% Allocate array sizes by n steps
x1=zeros(n);
x2=zeros(n);
V=zeros(n);
Z=zeros(n);

% Initial Conditions
x1(1)=-7;
x2(1)=7;
V(1)=-3;
Z(1)=3;

% Create differential equations
dVdt=@(x1,x2,V,Z) 1/m1*(-(k1+k2)*x1+k2*x2-(c1+c2)*V+c2*Z);
dZdt=@(x1,x2,V,Z) 1/m2*(k2*x1-(k2+k3)*x2+c2*V-(c2+c3)*Z);



for i=1:N-1
    V(i+1)=V(i)+h*dVdt(x1(i),x2(i),V(i),Z(i));
    x1(i+1)=x1(i)+h*V(i);
    Z(i+1)=Z(i)+h*dZdt(x1(i),x2(i),V(i),Z(i));
    x2(i+1)=x2(i)+h*Z(i);

end
figure(1)
plot(t,x1, 'blue',t,x2,'red');
xlabel('t');
ylabel('x');

fun = @(t,y)[y(2);(-(k1+k2)*y(1)+k2*y(3)-(c1+c2)*y(2)+c2*y(4))/m1;y(4);(k2*y(1)-(k2+k3)*y(3)+c2*y(2)-(c2+c3)*y(4))/m2];
tspan = [0 100];
y0 = [-7 -3 7 3];
[T ,Y] = ode45(fun,t,y0,odeset('RelTol',1e-12,'AbsTol',1e-12));
figure(2)
plot(T,[Y(:,1),Y(:,3)])

figure(3)
plot(t,[Y(:,1).'-x1;Y(:,3).'-x2])