Using Euler method to solve second order ODE

Hello,

I'm trying to write a program that usese Euler method to solve second order ODE. (represnts the movment of a package hanging from the roof with a spring).

I eventually will need to solve a more complex one, but as long as the simple one doesn't work, I have no reason to continue.

My ODE is: y'' = -25*y + 440.19 ;

For a reason I don't understand, my plot is not converging with the calculated solution for the ODE.

The code and the plot I'm getting are in the pictures. (the expected result was done with Runge-Kutta but since I can't manege to use it so solve a second order ODE I'm trying to understand the Euler method).

Zp_dotaim = y'' ; Zp_dot = y' ; Zp = y

Initial values : y'(0) = 0 ; y(0) = 18 ;

Thank you!

Zp_dot = zeros(1,1000) ; % prealocating vectors [m/sec]

Zp = zeros(1,1000) ; % prealocating vectors [m]

Zp(1) = 18 ; % initial condition

t = linspace(1,100,1000) ; % creating time vector [sec]

h = 0.01 ; % time step

for i = 1:999

% Zp_dotaim =@(Zp) 25*((2/(20-Zp))-1)*(Zp-20)-9.81 ;

Zp_dotaim =@(Zp) -25*Zp + 440.19 ;

Zp_dot(i+1) = Zp_dot(i) + h*Zp_dotaim(Zp(i)) ;

Zp(i+1) = Zp(i) + h*Zp_dot(i) ;

end

exact = (0.3924*cos(5*t)+17.6076);

plot(t,Zp,'r',t,exact,'b--') ; title('Zp(t)') ; grid on ;

xlabel('t [sec]') ; ylabel('Zp(t) [m]') ;

Answer 1

[Sam Chak]

Hi @Yuval Levy

Some minor issues. If you use Euler's method and you want to simulate for a relatively long duration, then you need to make the time step smaller. Also, the time vector begins from , because the initial value from the exact solution (cosine) is .

Zp_dot = zeros(1, 100001);      % prealocating vectors [m/sec]
Zp     = zeros(1, 100001);      % prealocating vectors [m]
Zp(1)  = 18;                    % initial condition
t = linspace(0, 10, 100001);    % creating time vector [sec]
h = 1e-4;                       % time step

for i = 1:100000
   Zp_dotaim   = @(Zp) - 25*Zp + 440.19 ;
   Zp_dot(i+1) = Zp_dot(i) + h*Zp_dotaim(Zp(i)) ; 
   Zp(i+1)     = Zp(i) + h*Zp_dot(i) ;
end

exact = (0.3924*cos(5*t)+17.6076);

plot(t, Zp, 'r', t, exact,'b--'); 
title('Zp(t)'); 
grid on;
xlabel('t [sec]');
ylabel('Zp(t) [m]') ; 

[Yuval Levy]

If you only knew how much time I spent trying to figure out what is wrong...

Thank you so much!

[Yuval Levy]

@Sam Chak

Follow up question- how did you choose the specific values for h and vector lengths? Is there a mathematical formula to do so?

Because it seems that every change I'm making ruins the results and I can't figure out the rule...

To be more specific- I need the final time to be 100 seconds, and not 10.

Thank you.

[Torsten]

If you work with constant stepsize, you can only pray that the chosen h is ok. Usually, you choose a value, calculate, half the value and compare the results. As soon as the results are no longer changing, you can respire.

Zp_dot = zeros(1, 10000001);      % prealocating vectors [m/sec]
Zp     = zeros(1, 10000001);      % prealocating vectors [m]
Zp(1)  = 18;                    % initial condition
t = linspace(0, 100, 10000001);    % creating time vector [sec]
h = 1e-5;                       % time step
for i = 1:10000000
    Zp_dotaim   = @(Zp) - 25*Zp + 440.19 ;
    Zp_dot(i+1) = Zp_dot(i) + h*Zp_dotaim(Zp(i)) ; 
    Zp(i+1)     = Zp(i) + h*Zp_dot(i) ;
end
exact = (0.3924*cos(5*t)+17.6076);
plot(t, Zp, 'r', t, exact,'b--'); 
title('Zp(t)'); 
grid on;
xlabel('t [sec]');
ylabel('Zp(t) [m]') ; 

[Sam Chak]

@Yuval Levy, I chose a smaller step size because I understand the limits of Euler's method.

The error in Euler's method tends to accumulate like a Snowball Rolling Downhill when the step size is relatively large in the fast-changing dynamics. This happens because the Euler's method uses the linear approximation approach on the dynamics, trucating higher-order terms.

In this example, I can use a relatively large step size because the dynamics is rather slow.

If you find my explanation is helpful, please consider accepting ✔ and voting 👍 the Answer. Thanks!

Zp_dot = zeros(1, 6001);      % prealocating vectors [m/sec]
Zp     = zeros(1, 6001);      % prealocating vectors [m]
Zp(1)  = 1;                    % initial condition
t = linspace(0, 60, 6001);    % creating time vector [sec]
h = 1e-2;                       % time step
for i = 1:6000
   Zp_dotaim   = @(Zp) - 0.01*Zp;
   Zp_dot(i+1) = Zp_dot(i) + h*Zp_dotaim(Zp(i)) ; 
   Zp(i+1)     = Zp(i) + h*Zp_dot(i) ;
end
exact = 1*cos(0.1*t);
plot(t, Zp, 'r', t, exact,'b--'); 
title('Zp(t)'); 
grid on;
xlabel('t [sec]');
ylabel('Zp(t) [m]') ; 

[Yuval Levy]

@Torsten Thank's alot! it worked. still don't really know how to figure out the h values and vector sizes, but you answed my question.

Thank you again.

[Torsten]

@Yuval Levy

still don't really know how to figure out the h values and vector sizes, but you answed my question.

I answered your question about the h value :-).

And the sizes of the vectors must be 1 x [(tend-tstart)/h + 1] in order that the solution and derivatives for all requested times can be saved within them.

Answer 2

[Tom Brenner]

You can't solve a second order differential equation with a single initial condition. You must have two. In this case, you assume that the first value of Zp_dot is zero (the vector was initialized with zeros) and add to this zero the approximate value of h times the second derivative.

You should determine how the exact solution 0.3924*cos(5*t)+17.6076 was arrived at (i.e., what the two initial conditions should be), and then correct your code.

[Yuval Levy]

Thanks for you answer Tom.

You are right, my bad. The initial value for Zp_dot is actually zero, I just forgot to write it in the explenation.

So unfortunately, that's not the origin of my problem... I have edited the post.