odesimp.m

很多matlab的源代码

function [tout, yout] = odesimp(F, t0, tf, y0, h)
% ODESIMP Numerical integration of state equations using Simpson's rule.
%
%	[T,Y] = ODESIMP('F',T0,TF,Y0,H) Numerical integration of state eqations 
%	using Simpon's Rule with fixed step size. 
%
%         T  - time points (column vector)
%         Y  - returned solution for states (column vectors)
%         F  - string function of t and y or  a function m-file F(t,y)
%         T0 - initial time
%         TF - final time
%         Y0 - initial condition vector
%         H  - step size. [Default: (TF-T0)/100]
%
%         ODESIMP (with no arguments) invokes the following example:
%
%         % Find and plot response of state equations in ODEMFILE.M for 0 to 2s 
%         % using step size = 0.05  and initial conditions Y0 = [3,-1]
%         >>[t1,y1] = odesimp('odemfile',0,2,[-3;1],0.05);
%         >>plot(t1,y1)


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help odesimp,disp('Strike a key to see results of the example')
pause,[t1,y1]=odesimp('odemfile',0,2,[-3;1],0.05);
vx=matverch;
if vx < 4, eval('clg');else,eval('clf');end
plot(t1,y1),return,end
 
vx=matverch;
flag=0;
if nargin < 5, h=(tf-t0)/100; end
t=t0;y=y0(:).';tout=t;yout=y;[m,n]=size(F);
if m==1 & n<9
if exist(F)>1,flag=1;
if vx<5
if exist('odesf')==3,
yout=eval('odesf(F,t0,tf,y0,h)');
[mm,nn]=size(yout);
tout=t0+(0:mm-1)*h;tout=tout(:);
return,
end
end
end
end
if flag==1;s=[F '(t,y)'];end

while t < tf
ta=t;ya=y;
if flag==1,s1=eval(s);s1=s1(:).';else,for j=1:m,s1(j)=eval(F(j,:));end,end
t=ta+h;y=ya+h*s1;
if flag==1,s2=eval(s);s2=s2(:).';else,for j=1:m,s2(j)=eval(F(j,:));end,end
t=ta+h/2;y=ya+h*(s1+s2)/4;
if flag==1,s3=eval(s);s3=s3(:).';else,for j=1:m,s3(j)=eval(F(j,:));end,end
t=ta;y=ya;t=t+h;y=y+h*(s1+4*s3+s2)/6;tout=[tout;t];yout=[yout;y];
end