matlab ode45 函数传自定义参数用法及定步长ode45结算函数

原文地址：matlab ode45 函数传自定义参数用法及定步长ode5解算函数 作者：jlxiaohuo

要用的时候总是忘记，这回给把它写在这里！

 

%%程序1

arg1 = 2;
arg2 = 1;
[T,Y] = ode45('vdp1000',[0 10],[2 0], [], arg1, arg2);
plot(T,Y(:,1),'-o');

 

%%程序2

function dy = vdp1000(t, y, flag, arg1, arg2)
dy =zeros(2,1);    %a column vector
dy(1) = y(2);
dy(2) = arg1*(arg2 - y(1)^2)*y(2) - y(1);

 

%%ode5

function Y = ode5(odefun,tspan,y0,varargin)

%ODE5 Solve differential equations with a non-adaptive method oforder 5.

% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN]integrates
% the system of differential equations y' = f(t,y) by stepping fromT0 to
% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a columnvector.
% The vector Y0 is the initial conditions at T0. Each row in thesolution
% array Y corresponds to a time specified in TSPAN.
%
% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additionalparameters
% P1,P2... to the derivative function asODEFUN(T,Y,P1,P2...).
% This is a non-adaptive solver. The step sequence is determined byTSPAN
% but the derivative function ODEFUN is evaluated multiple timesper step.
% The solver implements the Dormand-Prince method of order 5 in ageneral
% framework of explicit Runge-Kutta methods.
%
% Example
% tspan = 0:0.1:20;
% y = ode5(@vdp1,tspan,[2 0]);
% plot(tspan,y(:,1));
% solves the system y' = vdp1(t,y) with a constant step size of0.1,
% and plots the first component of the solution.

if ~isnumeric(tspan)
   
    error('TSPANshould be a vector of integration steps.');
   
end

if ~isnumeric(y0)
   
    error('Y0should be a vector of initial conditions.');
   
end

h = diff(tspan);

if any(sign(h(1))*h <= 0)
   
   error('Entries of TSPAN are not in order.')
   
end

try
   
    f0 =feval_r(odefun,tspan(1),y0,varargin{:});
   
catch
   
    msg =['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
   
   error(msg);
   
end

y0 = y0(:); % Make a column vector.

if ~isequal(size(y0),size(f0))
   
   error('Inconsistent sizes of Y0 and f(t0,y0).');
   
end

neq = length(y0);

N = length(tspan);

Y = zeros(neq,N);

% Method coefficients -- Butcher's tableau

%

% C | A

% --+---

% | B

C = [1/5; 3/10; 4/5; 8/9; 1];

A = [ 1/5, 0, 0, 0, 0
   
3/40, 9/40, 0, 0, 0

44/45 -56/15, 32/9, 0, 0

19372/6561, -25360/2187, 64448/6561, -212/729, 0

9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];

B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];

% More convenient storage

A = A.';

B = B(:);

nstages = length(B);

F = zeros(neq,nstages);

Y(:,1) = y0;

for i = 2:N
   
    ti =tspan(i-1);
   
    hi =h(i-1);
   
    yi =Y(:,i-1);
   
    % Generalexplicit Runge-Kutta framework
   
    F(:,1) =feval_r(odefun,ti,yi,varargin{:});
   
    for stage =2:nstages
       
       tstage = ti + C(stage-1)*hi;
       
       ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));
       
       F(:,stage) = feval_r(odefun,tstage,ystage,varargin{:});
       
    end
   
    Y(:,i) = yi+ F*(hi*B);
   
end

Y = Y.';