请教关于 Ode5
<P>定步长积分，听说可用Ode5来解，Ode5的用法和Ode45一样吗，请问哪里有Ode5下，谢谢大家了。</P>[此贴子已经被作者于2006-5-10 16:35:49编辑过]
回复：（qilin70007）请教关于 Ode5
<P>function Y = ode5(odefun,tspan,y0,varargin)<br>%ODE5Solve differential equations with a non-adaptive method of order 5.<br>% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = 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,);<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 <= 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 = ;</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 = ;</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编辑过]
谢谢
页:
[1]