Using Euler method to solve second order ODE

Question

Yuval Levy 2022-7-19

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1763625-using-euler-method-to-solve-second-order-ode

编辑： Torsten 2022-7-20

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]') ;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2022-7-19

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1763625-using-euler-method-to-solve-second-order-ode#answer_1010940

编辑：Sam Chak 2022-7-19

在 MATLAB Online 中打开

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]') ;

6 个评论
显示 4更早的评论隐藏 4更早的评论

Sam Chak 2022-7-20

编辑：Sam Chak 2022-7-20

在 MATLAB Online 中打开

@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 2022-7-20

@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 2022-7-20

编辑：Torsten 2022-7-20

@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 2022-7-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1763625-using-euler-method-to-solve-second-order-ode#answer_1010915

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Yuval Levy 2022-7-19

编辑：Yuval Levy 2022-7-19

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.

请先登录，再进行评论。

Using Euler method to solve second order ODE

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Using Euler method to solve second order ODE

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论