matlab变步长定步长_matlab optimvar「建议收藏」

概要

本文介绍matlab内置函数ode45中所采用的变步长策略以及相应部分代码的解读

笔者所使用的matlab版本为MTALAB R2018b

参考资料

matlab常微分方程求解器技术文档

The MATLAB ODE Suite, L. F. Shampine and M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997

DP54格式原论文

J. R. DORMAND AND P. J. PRINCE, A family of embedded Runge-Kutta formulae, J. Comp. Appl. Math., 6 (1980), pp. 19-26.

代码解读

ode45函数头

% ode为方程右端项，tspan为求解区间，y0为初值，options为odeset返回的参数结构体
function varargout = ode45(ode,tspan,y0,options,varargin)

初始参数设置

solver_name = 'ode45';  % 求解器名称
% Stats
nsteps  = 0;  % 求解步数
nfailed = 0;  % 重新迭代次数
nfevals = 0;  % 函数估值次数
% Output
FcnHandlesUsed  = isa(ode,'function_handle');  % 是否使用函数句柄
output_sol = (FcnHandlesUsed && (nargout==1));%sol=odeXX(...) 返回函数句柄
sol=[];
output_ty  = (~output_sol && (nargout > 0));%[t,y,...]=odeXX(...) 返回参数矩阵
% There might be no output requested...
% Handle solver arguments
% 参数设置
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
  options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, dataType] = ...
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

转至odearguments函数：

function [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, args, odeFcn, ...
          options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, ...
          dataType ] =   ...
    odearguments(FcnHandlesUsed, solver, ode, tspan, y0, options, extras)

  tspan = tspan(:);
  ntspan = length(tspan);
  if ntspan == 1    % Integrate from 0 to tspan 
    t0 = 0;          
    next = 1;       % Next entry in tspan.
  else              
    t0 = tspan(1);  
    next = 2;       % next entry in tspan
  end
  htspan = abs(tspan(next) - t0);
  tfinal = tspan(end);   
  y0 = y0(:);
  neq = length(y0);
  tdir = sign(tfinal - t0);
  f0 = feval(ode,t0,y0,args{ 
   :});
  dataType = superiorfloat(t0,y0,f0);  % 最高精度类型
  % Get the error control options, and set defaults.
  rtol = odeget(options,'RelTol',1e-3,'fast');  % 相对误差设置
  if rtol < 100 * eps(dataType) 
     rtol = 100 * eps(dataType);  % 误差不能小于100倍的当前类型可识别最小精度
  end
  atol = odeget(options,'AbsTol',1e-6,'fast');  % 绝对误差设置
  atol = atol(:);
threshold = atol / rtol;  % 绝对误差和相对误差比值，方便计算使用
% By default, hmax is 1/10 of the interval.
% 最大步长设置，默认是间隔的1/10
safehmax = 16.0*eps(dataType)*max(abs(t0),abs(tfinal));  % 数值安全步长
defaulthmax = max(0.1*(abs(tfinal-t0)), safehmax);
hmax = min(abs(tfinal-t0), abs(odeget(options,'MaxStep',defaulthmax,'fast')));
htry = abs(odeget(options,'InitialStep',[],'fast'));  % 初始步长设置
odeFcn = ode;
end

回到ode45函数：

refine = max(1,odeget(options,'Refine',4,'fast'));  % 插值个数 = refine-1
if ntspan > 2
  outputAt = 1;          % output only at tspan points 只在指定点需要插值
elseif refine <= 1
  outputAt = 2;          % computed points, no refinement 不需要插值
else
  outputAt = 3;          % computed points, with refinement 需要插值
  S = (1:refine-1) / refine;
end
t = t0;
y = y0;
nout = 0;
tout = []; yout = [];  % 输出参数初始化
if nargout > 0
    if ntspan > 2    % output only at tspan points 只在指定点输出值
      tout = zeros(1,ntspan,dataType);
      yout = zeros(neq,ntspan,dataType);
    else                                  % alloc in chunks
      chunk = min(max(100,50*refine), refine+floor((2^13)/neq));
      tout = zeros(1,chunk,dataType);
      yout = zeros(neq,chunk,dataType);
    end
  nout = 1;
  tout(nout) = t;
  yout(:,nout) = y;
end

注意上述代码中的refine参数，matlab的输出并不是全是由相应格式前进一步得到的，有一部分的解（默认为3/4）是通过插值公式得到的

DP54格式参数设置

% Initialize method parameters.
pow = 1/5;
A = [1/5, 3/10, 4/5, 8/9, 1, 1]; % Still used by restarting criteria
% B = [
% 1/5 3/40 44/45 19372/6561 9017/3168 35/384
% 0 9/40 -56/15 -25360/2187 -355/33 0
% 0 0 32/9 64448/6561 46732/5247 500/1113
% 0 0 0 -212/729 49/176 125/192
% 0 0 0 0 -5103/18656 -2187/6784
% 0 0 0 0 0 11/84
% 0 0 0 0 0 0
% ];
% E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];

% Same values as above extracted as scalars (1 and 0 values ommitted)
a2=cast(1/5,dataType);  % 转换数据类型
a3=cast(3/10,dataType);
a4=cast(4/5,dataType);
a5=cast(8/9,dataType);

b11=cast(1/5,dataType); 
b21=cast(3/40,dataType); 
b31=cast(44/45,dataType);
b41=cast(19372/6561,dataType);
b51=cast(9017/3168,dataType);
b61=cast(35/384,dataType);
b22=cast(9/40,dataType);
b32=cast(-56/15,dataType);
b42=cast(-25360/2187,dataType);
b52=cast(-355/33,dataType);
b33=cast(32/9,dataType);
b43=cast(64448/6561,dataType);
b53=cast(46732/5247,dataType);
b63=cast(500/1113,dataType);
b44=cast(-212/729,dataType);
b54=cast(49/176,dataType);
b64=cast(125/192,dataType);
b55=cast(-5103/18656,dataType);
b65=cast(-2187/6784,dataType);
b66=cast(11/84,dataType);

e1=cast(71/57600,dataType);
e3=cast(-71/16695,dataType);
e4=cast(71/1920,dataType);
e5=cast(-17253/339200,dataType);
e6=cast(22/525,dataType);
e7=cast(-1/40,dataType);

初始步长设置

hmin = 16*eps(t);  % 最小步长
if isempty(htry)
  % Compute an initial step size h using y'(t).
  % 根据初始点的右端项的值计算初始步长
  absh = min(hmax, htspan);
  rh = norm(f0 ./ max(abs(y),threshold),inf) / (0.8 * rtol^pow);
  if absh * rh > 1
    absh = 1 / rh;
  end
  absh = max(absh, hmin);
else
  absh = min(hmax, max(hmin, htry));
end
f1 = f0;

上述初始步长一般为$$0.8\ast \frac{max(绝对误差,相对误差\ast y)}{f\ast {相对误差}^{\frac{4}{5}}}$$

DP54格式计算

% Cleanup the main ode function call 
FcnUsed = isa(odeFcn,'function_handle');
odeFcn_main = odefcncleanup(FcnUsed,odeFcn,{ 
   });
% THE MAIN LOOP

done = false;
while ~done
% By default, hmin is a small number such that t+hmin is only slightly
  % different than t. It might be 0 if t is 0.
  hmin = 16*eps(t);
  absh = min(hmax, max(hmin, absh));    % couldn't limit absh until new hmin
  h = tdir * absh;   % 当前步长
  % Stretch the step if within 10% of tfinal-t.
  %是否抵达求解区间右端点
  if 1.1*absh >= abs(tfinal - t)
    h = tfinal - t;
    absh = abs(h);
    done = true;
  end
   % LOOP FOR ADVANCING ONE STEP.
   %单步循环
  nofailed = true;                      % no failed attempts
  while true
    y2 = y + h .* (b11.*f1 );
    t2 = t + h .* a2;
    f2 = odeFcn_main(t2, y2);
        
    y3 = y + h .* (b21.*f1 + b22.*f2 );
    t3 = t + h .* a3;
    f3 = odeFcn_main(t3, y3);
        
    y4 = y + h .* (b31.*f1 + b32.*f2 + b33.*f3 );
    t4 = t + h .* a4;
    f4 = odeFcn_main(t4, y4);
        
    y5 = y + h .* (b41.*f1 + b42.*f2 + b43.*f3 + b44.*f4 );
    t5 = t + h .* a5;
    f5 = odeFcn_main(t5, y5);
       
    y6 = y + h .* (b51.*f1 + b52.*f2 + b53.*f3 + b54.*f4 + b55.*f5 );
    t6 = t + h;
    f6 = odeFcn_main(t6, y6);

    tnew = t + h;
    if done
      tnew = tfinal;   % Hit end point exactly.
    end
    h = tnew - t;      % Purify h. 
    
    ynew = y + h.* ( b61.*f1 + b63.*f3 + b64.*f4 + b65.*f5 + b66.*f6 ); % 计算所得的数值解
    f7 = odeFcn_main(tnew,ynew);    
    
    nfevals = nfevals + 6;     % 进行了6次估值 

误差估计

 % Estimate the error.
    NNrejectStep = false;
    fE = f1*e1 + f3*e3 + f4*e4 + f5*e5 + f6*e6 + f7*e7;  % 误差估计
    err = absh * norm((fE) ./ max(max(abs(y),abs(ynew)),threshold),inf);                
     % Accept the solution only if the weighted error is no more than the
    % tolerance rtol. Estimate an h that will yield an error of rtol on
    % the next step or the next try at taking this step, as the case may be,
    % and use 0.8 of this value to avoid failures.
    if err > rtol                       % Failed step
      nfailed = nfailed + 1;  

不接受此步长的条件为$$h\ast f_E>max(绝对误差,相对误差\ast max(y_n,y_{n+1}))$$

更新步长

      if nofailed
         nofailed = false;
         absh = max(hmin, absh * max(0.1, 0.8*(rtol/err)^pow));
      else
        absh = max(hmin, 0.5 * absh);
      end
      h = tdir * absh;
      done = false;

步长更新公式为$$h_{new}=h\ast 0.8\ast (\frac{max(绝对误差,相对误差\ast max(y_n,y_{n+1}))}{h\ast f_E}{)}^{\frac{1}{5}}$$
最小更新为$$h_{new}=h\ast 0.1$$
若步长仍不被接受，则$$h_{new}=h\ast 0.5$$

接受步长

else                               % Successful step 
      break;   % 跳出while true循环
    end 
  end
nsteps = nsteps + 1;
% 生成输出
if output_ty || haveOutputFcn   
    switch outputAt
     case 2      % computed points, no refinement 不进行插值，直接输出
      nout_new = 1;
      tout_new = tnew;
      yout_new = ynew;
     case 3      % computed points, with refinement 进行插值
      tref = t + (tnew-t)*S;
      nout_new = refine;
      tout_new = [tref, tnew];
      yntrp45 = ntrp45split(tref,t,y,h,f1,f3,f4,f5,f6,f7,[]);  % ode45插值函数
      yout_new = [yntrp45, ynew];
     case 1      % output only at tspan points 只在特定点进行插值输出
      nout_new =  0;
      tout_new = [];
      yout_new = [];
      while next <= ntspan 
      if tdir * (tnew - tspan(next)) < 0
          break;                    % 当求值点越过输出点时进行插值
      end
        nout_new = nout_new + 1;              
        tout_new = [tout_new, tspan(next)];
        if tspan(next) == tnew
          yout_new = [yout_new, ynew];            
        else
          yntrp45 = ntrp45split(tspan(next),t,y,h,f1,f3,f4,f5,f6,f7,[]);
          yout_new = [yout_new, yntrp45];
        end  
        next = next + 1;
      end
    end
     % 写入输出
    if nout_new > 0   
      if output_ty
        oldnout = nout;
        nout = nout + nout_new;
        if nout > length(tout)
          tout = [tout, zeros(1,chunk,dataType)];  % requires chunk >= refine
          yout = [yout, zeros(neq,chunk,dataType)];
        end
        idx = (oldnout+1):nout;        
        tout(idx) = tout_new;
        yout(:,idx) = yout_new;
      end
    end  
  end 
  % If there were no failures compute a new h.
  % 当此步长被接受时仍需更改步长
  if nofailed   
    % Note that absh may shrink by 0.8, and that err may be 0.
    temp = 1.25*(err/rtol)^pow;
    if temp > 0.2
      absh = absh / temp;
    else
      absh = 5.0*absh;
    end
  end

此处的步长更新公式同上，但要求步长最多增大五倍

结束单步循环

% Advance the integration one step.
  t = tnew;
  y = ynew;
  f1 = f7;  % Already have f(tnew,ynew) 
end

程序结束

solver_output = { 
   };
if (nout > 0) % produce output
  if isempty(sol) % output [t,y,...]
    solver_output{ 
   1} = tout(1:nout).';
    solver_output{ 
   2} = yout(:,1:nout).';
    solver_output{ 
   end+1} = [nsteps, nfailed, nfevals];  % Column vector
   end
end
if nargout > 0
  varargout = solver_output;  % ode45输出
end 

细节回顾

• 必要参数有方程右端项(ode)，初始值(y0)，求解区间(tspan)
• 可选参数有绝对误差(atol,默认1e-6)，相对误差(rtol,默认1e-3)，最大步长(hmax,默认为间隔的1/10)，初始步长(htry)，插值标识(refine,默认为4)
• 当未设置初始步长时会根据初始右端项的值计算一个初始步长
• 会根据插值标识选择是否进行插值
• 步长更新倍率限制在0.1到5.0之间
• 步长更新公式$$h_{new}=h\ast 0.8\ast (\frac{max(绝对误差,相对误差\ast max(y_n,y_{n+1}))}{h\ast f_E}{)}^{\frac{1}{5}}$$

小结

本文只保留了ode45源代码中的变步长及一般常微分方程求解部分，去除了其他无关变步长的内容（比如其他可选参数的设置及其作用；形如$M(t)y^{\prime }=f(x,y)$常微分方程的求解等等），但这并不影响上述代码变步长部分的运行。ode45的变步长策略并不复杂，参考资料中的技术文档发行年份是1997年，其中的大部分篇幅是在讲matlab在求解刚性问题中所做的努力，ode45在其中只占很小一部分，甚至变步长策略都没有提及，DP54格式原论文的发表年份是1980年，文中提到的变步长策略和现在使用的已经十分相似。ode45在求解方程时会通过插值来提高求解效率和绘图的光滑性，对于求解精度并没有太过严格的要求。

补充

上述代码中出现了两个matlab的私有函数（odefcncleanup和ntrp45split），用户没法直接调用，可以将它们拷贝到用户目录下来进行调用。

今天的文章matlab变步长定步长_matlab optimvar「建议收藏」分享到此就结束了，感谢您的阅读。