Fourth Order Runge Kutta for Systems

9 views (last 30 days)
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];
% 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;
% 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);

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

More Answers (2)

James Tursa
James Tursa on 18 Apr 2019
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;
James Tursa
James Tursa on 18 Apr 2019
Edited: 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;
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.

Sign in to comment.

Meysam Mahooti
Meysam Mahooti on 5 May 2021


Find more on Startup and Shutdown in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!