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]') ;
采纳的回答
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 
.
, 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]') ;
6 个评论
Yuval Levy
2022-7-19
Yuval Levy
2022-7-19
编辑：Yuval Levy
2022-7-19
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]') ;
@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.
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
2022-7-20
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.
更多回答（1 个）
Tom Brenner
2022-7-19
1 个投票
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Numeric Solvers 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
