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 · 2022-11-27T21:05:19.987

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

Revision 1

The error returned is the actual step error. The error that is required for the step size control however is the unit step error or error density, which is the step error with divided by h

  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/h;
  end

Changing the example to a simple bi-stable model oscillating between two branches of stable equilibria

  rhs = @(x,y) 3*y-y^3 + 3*cos(x);
  x0 = 0; y0 = 1; x_end = 13.5; tol = 5e-3; h_init = 5e-2;

gives plots of solution, error (against an ode45 integration) and step sizes

Red crosses are the step sizes of rejected steps.

Revision 2

The error in the function values can be used as an error guidance for the extrapolation value which is of 5th order, making the method a 5th order method in extrapolation mode. As it uses the 4th order error to predict the 5th order optimal step size, a caution factor is recommended, the code changes in the appropriate places to

    factor = 0.75*abs(tol/err)^0.2;


...

  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+(y1-y2)/15;
    err = (y1-y2)/15;
  end

In the plots the step size is appropriately larger, but the error shows sharper and larger spikes, this version of the method is apparently less stable.

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

1 Answers1

Revision 1

Revision 2

Linked