请教关于 Ode5

qilin70007

<P>定步长积分，听说可用Ode5来解，Ode5的用法和Ode45一样吗，请问哪里有Ode5下，谢谢大家了。</P>

[此贴子已经被作者于2006-5-10 16:35:49编辑过]

happy

<P>function Y = ode5(odefun,tspan,y0,varargin)<br>%ODE5  Solve differential equations with a non-adaptive method of order 5.<br>%   Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates <br>%   the system of differential equations y' = f(t,y) by stepping from T0 to <br>%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.<br>%   The vector Y0 is the initial conditions at T0. Each row in the solution <br>%   array Y corresponds to a time specified in TSPAN.<br>%<br>%   Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters <br>%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). <br>%<br>%   This is a non-adaptive solver. The step sequence is determined by TSPAN<br>%   but the derivative function ODEFUN is evaluated multiple times per step.<br>%   The solver implements the Dormand-Prince method of order 5 in a general <br>%   framework of explicit Runge-Kutta methods.<br>%<br>%   Example <br>%         tspan = 0:0.1:20;<br>%         y = ode5(@vdp1,tspan,[2 0]);  <br>%         plot(tspan,y(:,1));<br>%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, <br>%     and plots the first component of the solution.   </P>
<P>if ~isnumeric(tspan)<br>  error('TSPAN should be a vector of integration steps.');<br>end</P>
<P>if ~isnumeric(y0)<br>  error('Y0 should be a vector of initial conditions.');<br>end</P>
<P>h = diff(tspan);<br>if any(sign(h(1))*h &lt;= 0)<br>  error('Entries of TSPAN are not in order.') <br>end  </P>
<P>try<br>  f0 = feval(odefun,tspan(1),y0,varargin{:});<br>catch<br>  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];<br>  error(msg);  <br>end  </P>
<P>y0 = y0(:);   % Make a column vector.<br>if ~isequal(size(y0),size(f0))<br>  error('Inconsistent sizes of Y0 and f(t0,y0).');<br>end  </P>
<P>neq = length(y0);<br>N = length(tspan);<br>Y = zeros(neq,N);</P>
<P>% Method coefficients -- Butcher's tableau<br>%  <br>%   C | A<br>%   --+---<br>%     | B</P>
<P>C = [1/5; 3/10; 4/5; 8/9; 1];</P>
<P>A = [ 1/5,          0,           0,            0,         0<br>      3/40,         9/40,        0,            0,         0<br>      44/45        -56/15,       32/9,         0,         0<br>      19372/6561,  -25360/2187,  64448/6561,  -212/729,   0<br>      9017/3168,   -355/33,      46732/5247,   49/176,   -5103/18656];</P>
<P>B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];</P>
<P>% More convenient storage<br>A = A.'; <br>B = B(:);      </P>
<P>nstages = length(B);<br>F = zeros(neq,nstages);</P>
<P>Y(:,1) = y0;<br>for i = 2:N<br>  ti = tspan(i-1);<br>  hi = h(i-1);<br>  yi = Y(:,i-1);  <br>  <br>  % General explicit Runge-Kutta framework<br>  F(:,1) = feval(odefun,ti,yi,varargin{:});  <br>  for stage = 2:nstages<br>    tstage = ti + C(stage-1)*hi;<br>    ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));<br>    F(:,stage) = feval(odefun,tstage,ystage,varargin{:});<br>  end  <br>  Y(:,i) = yi + F*(hi*B);<br>  <br>end<br>Y = Y.';<br></P>

[此贴子已经被作者于2006-5-10 16:44:06编辑过]

qilin70007

谢谢[em01]