How to perform adaptive step size using Runge-Kutta fourth order (Matlab)?

Question

For me, it seems like the estimated hstep takes quite a long time and long iteration to converge. I tried it with this first ODE. Basically, you perform the difference between RK4 with stepsize of h with h/2.Please note that to reach the same timestep value, you will have to use the y value after two timestep of h/2 so that it reaches h also.

frhs=@(x,y) x.^2*y;

Is my code correct?

clear all;close all;clc
c=[]; i=1; U_saved=[]; y_array=[]; y_array_alt=[];
y_arr=1; y_arr_2=1;

frhs=@(x,y) 20*cos(x); 

tol=0.001;
y_ini= 1;
y_ini_2= 1;
c=abs(y_ini-y_ini_2)
hc=1
all_y_values=[];
for m=1:500
  if (c>tol || m==1)
      fprintf('More')
      y_arr
      [Unew]=vpa(Runge_Kutta(0,y_arr,frhs,hc))
      if (m>1)
         y_array(m)=vpa(Unew);
         y_array=y_array(logical(y_array));
         
      end
   
      [Unew_alt]=Runge_Kutta(0,y_arr_2,frhs,hc/2);
      [Unew_alt]=vpa(Runge_Kutta(hc/2,Unew_alt,frhs,hc/2))
      if (m>1)
          y_array_alt(m)=vpa(Unew_alt);
          y_array_alt=y_array_alt(logical(y_array_alt));
         
      end
      fprintf('More')
      %y_array_alt(m)=vpa(Unew_alt);
      c=vpa(abs(Unew_alt-Unew) )
      hc=abs(tol/c)^0.25*hc
      if (c<tol)
          fprintf('Less')
          
          y_arr=vpa(y_array(end) )
          y_arr_2=vpa(y_array_alt(end) )
          [Unew]=Runge_Kutta(0,y_arr,frhs,hc)
          all_y_values(m)=Unew;
          [Unew_alt]=Runge_Kutta(0,y_arr_2,frhs,hc/2);
          [Unew_alt]=Runge_Kutta(hc/2,Unew_alt,frhs,hc/2)
          c=vpa( abs(Unew_alt-Unew) )
          hc=abs(tol/c)^0.2*hc
      end
  end
end

all_y_values

Your right side `rhs` depends markedly on `x`, but you do not pass `x`, only the fixed values `0` and `h/2`. — Lutz Lehmann, Apr 12 '21 at 15:53

Lutz Lehmann · Answer 1 · 2021-04-14T08:16:50.833

A better structure for the time loop has only one place where the time step is computed.

  x_array = [x0]; y_array = [y0]; h = h_init;
  x = x0; y = y0;
  while x < x_end
    [y_new, err] = RK4_step_w_error(x,y,rhs,h);
    factor = abs(tol/err)^0.2;
    if factor >= 1
      y_array(end+1) = y = y_new;
      x_array(end+1) = x = x+h;
    end
    h = factor*h;
  end

For the data given in the code

  rhs = @(x,y) 20*cos(x);
  x0 = 0; y0 = 1; x_end = 6.5; tol = 1e-3; h_init = 1;

this gives the result against the exact solution

The computed points lie exactly on the exact solution, for the segments between them one would need to use a "dense output" interpolation. Or as a first improvement, just include the middle value from the half-step computation.

  function [ y_next, err] = RK4_step_w_error(x,y,rhs,h)
    y2 = RK4_step(x,y,rhs,h);
    y1 = RK4_step(x,y,rhs,h/2);
    y1 = RK4_step(x+h/2,y1,rhs,h/2);
    y_next = y1;
    err = (y2-y1)/15;
  end 
  
  function y_next = RK4_step(x,y,rhs,h)
    k1 = h*rhs(x,y);
    k2 = h*rhs(x+h/2,y+k1);
    k3 = h*rhs(x+h/2,y+k2);
    k4 = h*rhs(x+h,y+k3);
    y_next = y + (k1+2*k2+2*k3+k4)/6;
  end

How to perform adaptive step size using Runge-Kutta fourth order (Matlab)?

1 Answers1

Linked