3种微分方程求解器

Euler中点法、4阶龙格库塔法和对角隐式辛积分R-K方法，Matlab单步求解。

%%  p281页 对角隐式辛积分R-K方法 精度为4阶  1e-4
%  调用格式：[TZ] = DiagImpSym（functspanh0z0ers）
%  其中：    func为微分方程，i.e. func = @（tz）cos（t）
%            tspan为起止时间，tspan=[0 2*pi] h0为步长
%            z0为微分方程状态变量的初始值z0 = 0 ers为迭代精度误差上限
function [TZ] = DiagImpSym（functspanh0z0ers）
   %% 公式系数
%    b2 = -0.53652708;
%    b3 = 2.37893931;
%    b4 = 1.8606818856;
%    b1 = 1.0-b2-b3-b4;
   load DiagImpSymCoef.mat
   b1 = X（1）;
   b2 = X（2）;
   b3 = X（3）;
   b4 = X（4）;
   z0 = reshape（z0length（z0）1）;
   %% 时间节点
   t0 = tspan（1）;
   tF = tspan（end）;
   if（t0>=tF）
       error（‘起止时间输入错误！‘）;
   end
%    h0 = 1e-5;%步长
   if（h0>tF-t0）
       error（‘输入的步长过大！‘）
   end
   nh = floor（（tF-t0）/h0）;
   if（tF>=t0+nh*h0）
      T = linspace（t0tFnh+1）.‘;
      h = （tF-t0）/（nh+1）;
   else
      T = linspace（t0tFnh）.‘;
      h = （tF-t0）/nh;
   end
   num_T = length（T）;
  
   %% 微分方程的输入个数
   num_var = nargin（func）;
   if（num_var~=2）       
      error（‘微分方程的输入个数错误！‘）; 
   end
   %% 求解微分方程
   zk = z0;
   Z = zeros（num_Tlength（z0））;
   Z（1:） = zk;
   for k = 1:num_T-1   
       % t
       tk = T（k）;
       tk1 = T（k+1）;
       h = tk1-tk;
       % Y1
       Y10 = zk+b1*h/2*feval（functk+b1*h/2zk）;%时间节点的问题，初值问题
       while（1）
           Y1 = zk+b1*h/2*feval（functk+b1*h/2Y10）;
           if（norm（Y1-Y10）/norm（Y1）               break;
           else
               Y10 = Y1; 
           end           
       end
       % Y2
        Y20 = 2*Y1-zk+b2*h/2*feval（functk+（2*b1+b2）*h/22*Y1-zk）;%时间节点的问题，初值问题
       while（1）
           Y2 = 2*Y1-zk+b2*h/2*feval（functk+（2*b1+b2）*h/2Y20）;
           if（norm（Y2-Y20）/norm（Y2）               break;
           else
               Y20 = Y2; 
           end           
       end
       % Y3
        Y30 = 2*Y2-（2*Y1-zk）+b3*h/2*feval（functk+（2*b1+2*b2+b3）*h/22*Y2-（2*Y1-zk））;%时间节点的问题，初值问题
       while（1）
            Y3 = 2*Y2-（2*Y1-zk）+b3*h/2*feval（functk+（2*b1+2*b2+b3）*h/2Y30）;
           if（norm（Y3-Y30）/norm（Y3）               break;
           else
               Y30 = Y3; 
           end           
       end
       % Y4
        Y40 = 2*Y3-（2*Y2-2*Y1+zk）+b4*h/2*feval（functk+（2*b1+2*b2+2*b3+b4）*h/22*Y3-（2*Y2-2*Y1+zk））;%时间节点的问题，初值问题
       while（1）
            Y4 = 2*Y3-（2*Y2-2*Y1+zk）+b4*h/2*feval（functk+（2*b1+2*b2+2*b3+b4）*h/2Y40）;
           if（norm（Y4-Y40）/norm（Y4）               break;
           else
               Y40 = Y4; 
           end           
       end
       % zk1
       zk1 = 2*Y4-（2*Y3-2*Y2+2*Y1-zk）;
       Z（k+1:） = zk1;
       zk = zk1;
   end
end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2890  2015-07-16 21:45  三种求解器\DiagImpSym.m

     文件        208  2015-07-16 21:08  三种求解器\DiagImpSymCoef.mat

     文件       1225  2015-07-16 21:45  三种求解器\EulerCenterForm.m

     文件        902  2015-07-17 09:25  三种求解器\Example_1.m

     文件       1090  2015-07-16 21:45  三种求解器\R_K_4.m

     文件       2432  2015-07-16 21:09  三种求解器\SolveRoots.m

     目录          0  2015-07-17 09:47  三种求解器

----------- ---------  ---------- -----  ----

                 8747                    7