如何使用 Runge-Kutta 四阶(Matlab)执行自适应步长? [英] How to perform adaptive step size using Runge-Kutta fourth order (Matlab)?

问题描述

对我来说，估计的 hstep 似乎需要很长时间和长时间的迭代才能收敛.我用第一个 ODE 尝试过.基本上，您执行具有 h/2 步长的 RK4 之间的差异.请注意，要达到相同的时间步长值，您必须在 h/2 的两个时间步长后使用 y 值，以便它也达到 h.

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;

我的代码正确吗?

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

推荐答案

更好的时间循环结构只有一个地方可以计算时间步长.

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

对于代码中给出的数据

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

这给出了精确解的结果

计算点正好位于精确解上，对于它们之间的线段，需要使用密集输出"；插值.或者作为第一个改进，只包括半步计算的中间值.

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 

这篇关于如何使用 Runge-Kutta 四阶(Matlab)执行自适应步长?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！