gill.m

提供积分算法变步长Gill法代码

function [y_out t_out] = Gill(Fun,tspan,y_0,options,para)
% March 31, 2006
% Coded by Huabin XI.
% [e-mail]
% www.baisi.net

%==========================================================================
% This is the mathematical computational method of step-vary Gill method.
% 
% Instruction:
% All the varibles written in upper case are vector varibles
% Fun            : differential equations,m*1
% tspan          : time varible, [t_0 t_final]
% y_0            : initial values,m*1
% h              : step size
% epsilon        : absolute error tolerance
% L              : minimum step size, h/(2^L), integer varible
% options        : [h epsilon L]
% para           : struct varible, used for updating coefficients of
%                  differential equations
% y_out          : output
% t_out          : output, m*n, n: number of time pionts, m: number of states
%==========================================================================
% Check inputs
if nargin ==5 && isstruct(para)==1
                                          % Add user's define data  
else
    para    = []; 
end
%==========================================================================
t_0      = tspan(1);
                                          % begin time
t_final  = tspan(2);
                                          % final time                                       
h_0      = options(1);
                                          % intial step size                                          
epsilon  = options(2);
                                          % absolute error tolerance
L        = options(3);
                                          % stepsize limit
X        = 1;
                                          % counter                                      
q_0      = 0;
                                          % auxilary varible(used in the loop)
y_0      = y_0';

t        = t_0;   
                                          % (used in the loop)
h        = h_0;
                                          % (used in the loop)

y_out    = y_0;
t_out    = t_0;
%==========================================================================
% exit conditions :  reached minimum stepsize limit, success
while 1
    
    [y_a q_a]= V_G(Fun,h/2,t,y_0,q_0,para);
    [y_b q_b]= V_G(Fun,h/2,t+h/2,y_a,q_a,para);
                                              % two steps(h/2) integrations
    [y_c q_c]= V_G(Fun,h,t,y_0,q_0,para);
                                              % one step(h) interation
    if sum((y_b-y_c).^2)^0.5>epsilon
        h            = h/2;
        X            = 2*X;
                                              % descrease stepsize
        y_c          = y_b;
                                              % get two steps' value
        if X>=L
            disp('reach the minimum step size!');
            break;
                                              % reach minimum step size
        end
    else
        y_out        = [y_out y_b];
        t_out        = [t_out t+h];
        t            = t+h; 
        y_0          = y_b;
        q_0          = q_b;
    end
    if t>=t_final
        disp('differential equations have been solved!');
        break;
                                              %
    end
end

function [y_4 q_4]=V_G(Fun,h,t,y_0,q_0,para)
 % step 1
    K_1      = h*feval(Fun,t,y_0,para);
    y_1      = y_0+(K_1-2*q_0)/2;
    q_1      = q_0+3*(K_1-2*q_0)/2-K_1/2;
    % step 2
    K_2      = h*feval(Fun,t+h/2,y_1,para);
    y_2      = y_1+(1-(1/2)^0.5)*(K_2-q_1);
    q_2      = q_1+3*(1-(1/2)^0.5).*(K_2-q_1)-(1-(1/2)^0.5)*K_2;
    % step 3
    K_3      = h*feval(Fun,t+h/2,y_2,para);
    y_3      = y_2+(1+(1/2)^0.5)*(K_3-q_2);
    q_3      = q_2+3*(1+(1/2)^0.5)*(K_3-q_2)-(1+(1/2)^0.5)*K_3;
    % step 4
    K_4      = h*feval(Fun,t+h,y_3,para);
    y_4      = y_3+(K_4-2*q_3)/6;
    q_4      = q_3+3*(K_4-2*q_3)/6-K_4/2;