CODEWITHSUNDEEP
matlab
ak ve
ak ve

Reputation: 29

How to generalize the RK4 code to solve m number of ODEs in MATLAB?

The following MATLAB code represents the solution of two first order ODEs:


dy1 = @(x,y1,y2) y2; 
dy2 = @(x,y1,y2) 4*x + 10*sin(x) - y1;

%initial values
x0 = pi;
xn = 2*pi;
y1 = 0; % initial value of y1
y2 = 2; % initial value of y2  
h = 0.1; % step size

for x = x0 : h : xn-h
  L1 = h*dy2(x, y1, y2);
  K1 = h*dy1(x, y1, y2);
  L2 = h*dy2(x + h/2, y1 + K1/2, y2 + L1/2);
  K2 = h*dy1(x + h/2, y1 + K1/2, y2 + L1/2);
  L3 = h*dy2(x + h/2, y1 + K2/2, y2 + L2/2);
  K3 = h*dy1(x + h/2, y1 + K2/2, y2 + L2/2);
  L4 = h*dy2(x + h, y1 + K3, y2 + L3);
  K4 = h*dy1(x + h, y1 + K3, y2 + L3);
  y2 = y2 + (L1 + 2*L2 + 2*L3 + L4)/6;
  y1 = y1 + (K1 + 2*K2 + 2*K3 + K4)/6;
  x = x + h;
  fprintf ('%f \t %f\t %f\n',x,y2,y1);
end

How can I make this code general so it could easily manage to solve m number of decomposed first order ODEs (i.e., dy1, dy2, ....., dym)?

Many thanks for help.

Upvotes: 0

Views: 123

Answers (1)

James Tursa
James Tursa

Reputation: 2636

Just use state vectors as if you were coding for ode45( ) instead of individual scalar variables for your states. E.g.,

dy = @(x,y) [y(2); 4*x + 10*sin(x) - y(1)];
%initial values
x0 = pi;
xn = 2*pi;
y0 = [0;2]; % initial value of y
h = 0.1; % step size
x = x0 : h : xn-h;
n = numel(x);
y = zeros(numel(y0),n);
y(:,1) = y0;
for i=1:n-1
    K1 = h * dy(x(i)      , y(:,i)       );
    K2 = h * dy(x(i) + h/2, y(:,i) + K1/2);
    K3 = h * dy(x(i) + h/2, y(:,i) + K2/2);
    K4 = h * dy(x(i) + h  , y(:,i) + K3  );
    y(:,i+1) = y(:,i) + (K1 + 2*K2 + 2*K3 + K4)/6;
end

If you ever need to change the number of state variables and differential equations, you simply need to change the function handle and initial conditions. The looping code stays exactly the same. Also, you can pass the function handle, range, and initial conditions to ode45( ) and compare results directly with your results. E.g.,

[x,y] = ode45( dy, [x0 xn], y0 );

Upvotes: 1

Related Questions