MATLAB Answers

系统的第四阶runge kutta

48 views (last 30 days)
Show older comments
Angel Nunez
Angel Nunez on 18 Apr 2019
  • 0
Answered: Meysam Mahooti on 5 May 2021
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
百分比持有该系统的M文件。
功能dy = textbook_example_sys(t,y)
dy = [y(2); 9.81-0.25/68.1*y(2)^2];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% runge-kutta system code
功能[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;
fori = 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）;

Sign in to answer this question.

Accepted Answer

Angel Nunez
Angel Nunez on 18 Apr 2019
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)';

更多的Answers (2)

詹姆斯·图萨（James Tursa）
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;
应该是这个
y_vals（i+1，:) = y_vals（i，:)+（k_1+2*（k_2+k_3）+k_4）/6*h;
2 Comments
Show＆nbsp1旧评论
詹姆斯·图萨（James Tursa）
詹姆斯·图萨（James Tursa） on 18 Apr 2019
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;
我建议采用列矢量方法，使其与Ode45和朋友兼容。因此，将所有y_vals（i，:)等更改为y_vals（：：，i）等。

Sign in to comment.

Meysam Mahooti
Meysam Mahooti on 5 May 2021
//www.tatmou.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle
//www.tatmou.com/matlabcentral/fileexchange/61130-runge-kutta-fehlberg-rkf78?s_tid=srchtitle

Categories

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

开始狩猎！