Fourth Order Runge Kutta for Systems

Question

Angel Nunez 2019-4-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/457227-fourth-order-runge-kutta-for-systems

回答： Meysam Mahooti 2021-5-5

I wrote a functions that takes in a system of first order differential equations from a file named dydt_sys and solves them. I got the code to run but my numbers are off. Any help would be appreciated.. I should get

0 0 0

2.0000 19.1656 18.7256

4.0000 71.9311 33.0995

6.0000 147.9521 42.0547

8.0000 237.5104 46.9345

10.0000 334.1626 49.4027

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% seperate m-file that holds the system. 
function dy = textbook_example_sys(t,y)
dy = [y(2); 9.81-0.25/68.1*y(2)^2];
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% runge-kutta system code
function [t_vals,y_vals] = RK4_sys(dydt_sys,t_span,y_0,h)
 
t_start = t_span(1); 
t_final = t_span(2); 
t_vals = (t_start:h:t_final)';
num_steps = length(t_vals);
y_vals(1,:) = y_0; 
for i=1:num_steps-1  
    
    k_1 = dydt_sys(t_vals(i),y_vals(i,:)); 
    k_2 = dydt_sys(t_vals(i)+h/2,y_vals(i,:)+k_1/2*h); 
    k_3 = dydt_sys(t_vals(i)+h/2,y_vals(i,:)+k_2/2*h); 
    k_4 = dydt_sys(t_vals(i)+h,y_vals(i,:)+k_3*h); 
    
    y_vals(i+1,:) = y_vals(i) + (k_1+2*(k_2+k_3)+k_4)/6*h; 
    
end
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% command line
>> t_span = [0 10];y_0 = [0 0]; h = 2; 
>> [t_vals,y_vals] = RK4_sys(@textbook_example_sys,t_span,y_0,h); 

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

请先登录，再进行评论。

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

Answer 1

Angel Nunez 2019-4-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/457227-fourth-order-runge-kutta-for-systems#answer_371273

在 MATLAB Online 中打开

I still kept getting an error when i did that but it got me thinking about how I was playing with column and row vectors. I ended up transpoing all of them and it fixed the problem. Thank you!

    k_1 = dydt_sys(t_vals(i),y_vals(i,:))'; 
    k_2 = dydt_sys(t_vals(i)+h/2,y_vals(i,:)+k_1/2*h)'; 
    k_3 = dydt_sys(t_vals(i)+h/2,y_vals(i,:)+k_2/2*h)'; 
    k_4 = dydt_sys(t_vals(i)+h,y_vals(i,:)+k_3*h)'; 

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

请先登录，再进行评论。

Answer 2

James Tursa 2019-4-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/457227-fourth-order-runge-kutta-for-systems#answer_371248

在 MATLAB Online 中打开

This line does not have the rhs state syntax correct:

y_vals(i+1,:) = y_vals(i) + (k_1+2*(k_2+k_3)+k_4)/6*h;

It should be this instead

y_vals(i+1,:) = y_vals(i,:) + (k_1+2*(k_2+k_3)+k_4)/6*h;

2 个评论
显示无隐藏无

Angel Nunez 2019-4-18

That's what i originally thought, but I got the error

Unable to perform assignment because the size of the left side is 1-by-2 and the size of the right side is 2-by-2.

Error in RK4_sys (line 16)

y_vals(i+1,:) = y_vals(i,:) + (k_1+2*(k_2+k_3)+k_4)/6*h;

James Tursa 2019-4-18

编辑：James Tursa 2019-4-18

在 MATLAB Online 中打开

You need to decide whether your state vector is going to be a row vector or a column vector and be consistent. Currently you have it mixed. E.g. your derivative function is a column vector:

dy = [y(2); 9.81-0.25/68.1*y(2)^2];

but your state update equations are row vectors:

y_vals(i+1,:) = y_vals(i,:) + (k_1+2*(k_2+k_3)+k_4)/6*h;

I would suggest going with the column vector approach to make it compatible with ode45 and friends. So change all of the y_vals(i,:) etc to y_vals(:,i) etc.

请先登录，再进行评论。