MatLa求解延迟微分方程的注意事项.docx

1、MatLa求解延迟微分方程的注意事项MatLab求解延迟微分方程的注意事项 2010-08-14 17:04MatLab求解延迟微分方程的注意事项使用MATLAB求解延时微分方程的两种方法：DDE23和SimuLink有些不同点需要注意，否则结果会出现错误使用MATLAB来求解延迟微分方程是在生物数学和化学计算求解中经常遇到的事，在其它领域也比较常见。我所知道的，在MATLAB中求解微分方程有三种方法：1.使用ode45（龙格-库塔法的一个变种）求解，这时用一个数组，记录y的延迟项，但是c的值要考虑步长，再代入方程就能实现延时效应；2.使用dde23求解常数延时方程、使用ddesd可以用来求解

2、延迟与时间t有关的延迟微分方程；3.使用SimuLink建模求解，SimuLink是求解广义微分代数系统的通用工具，功能很强大，但是看惯了编程指令的人可能不大习惯，调试似乎也不太方便。 既然本文专门讨论求解延迟微分方程，就先介绍一下专用工具dde23，dde系列求解函数是由Southern Methodist University 的L.F. Shampine 和S. Thompson根据他们早期使用Fortran编写的Fortran 90 DDE Solver 移植到MATLAB上的，从MATLAB6.5开始加入MATLAB的官方发行版，根据薛定宇教授在其几本关于MATLAB的著作中提到的，

3、该函数返回的sol中结构体sol.x和sol.y均为按行排列，与ode45等不同，不太规范（没办法，因为这个函数本来就不是Mathworks的官方作品），不过这一点已经不大可能得到改进了，因为L.F. Shampine 和S. Thompson已经决定停止维护这个文件。如果您想进一步了解该函数，可以访问它的主页。 MATLAB帮助中关于该函数的介绍不很清楚，如果需要进一步了解这个函数，需要下载作者为其写的手册。下面以MATLAB中所附的一个例程来说明这个函数与Simulink建模求解的不同。 在MATLAB prompt中键入edit ddex1就会找看到函数作者所写的一个入门例子：funct

4、ion ddex1%DDEX1 Example 1 for DDE23.% This is a simple example of Wille and Baker that illustrates the% straightforward formulation, computation, and plotting of the solution % of a system of delay differential equations (DDEs). % The differential equations% y_1(t) = y_1(t-1) % y_2(t) = y_1(t-1)+y_2

5、(t-0.2)% y_3(t) = y_2(t)% are solved on 0, 5 with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for% t = 0. % The lags are specified as a vector 1, 0.2, the delay differential% equations are coded in the subfunction DDEX1DE, and the history is% evaluated by the function DDEX1HIST. Because the history i

6、s constant it% could be supplied as a vector:% sol = dde23(ddex1de,1, 0.2,ones(3,1),0, 5);% See also DDE23, FUNCTION_HANDLE.% Jacek Kierzenka, Lawrence F. Shampine and Skip Thompson% Copyright 1984-2004 The MathWorks, Inc.% $Revision: 1.2.4.2 $ $Date: 2005/06/21 19:24:16 $sol = dde23(ddex1de,1, 0.2,

7、ddex1hist,0, 5);figure;plot(sol.x,sol.y)title(An example of Wille and Baker.); xlabel(time t);ylabel(solution y);% -function s = ddex1hist(t)% Constant history function for DDEX1.s = ones(3,1);% -function dydt = ddex1de(t,y,Z)% Differential equations function for DDEX1.ylag1 = Z(:,1);ylag2 = Z(:,2);

8、dydt = ylag1(1)ylag1(1) + ylag2(2)y(2) ; 这里先不管函数使用的具体语法，求解模型为：显然有两个延时常数1、0.2。在使用手册中这样介绍：the DDE reduces to an initial value problem for an ODE with y(t 1) equal to the given history S(t 1) and initial value y(0) = 1也就是说，仿真求解开始时间为t=0，此时y_1=1，要特别注意，这里设定的y的history并不仅仅是y的初始值，也包括t tau=1;%给定延迟时间 history=0

9、;%初始值 tapan0,10;%求解时间范围 sol=dde23(yanchi,tau,history,tapan);%延迟问题求解 plot(sol.x,sol.y)；%作图下面附上了图片x(0)=0和x(0)=2的情况显然初始值不同结果不同，这就是为什么需要初始值的情况延迟微分方程是：dx(t)/dt=0.2*x(t-17)/(1+x(t-17)10)-0.1x(t)初始条件是x(0)=1.2当t0 then each ioutp-th point will be print. % % Output parameters: % Texp - time values % Lexp - Ly

10、apunov exponents to each time value. % % Users have to write their own ODE functions for their specified % systems and use handle of this function as rhs_ext_fcn - parameter. % % Example. Lorenz system: % dx/dt = sigma*(y - x) = f1 % dy/dt = r*x - y - x*z = f2 % dz/dt = x*y - b*z = f3 % % The Jacobi

11、an of system: % | -sigma sigma 0 | % J = | r-z -1 -x | % | y x -b | % % Then, the variational equation has a form: % % F = J*Y % where Y is a square matrix with the same dimension as J. % Corresponding m-file: % function f=lorenz_ext(t,X) % SIGMA = 10; R = 28; BETA = 8/3; % x=X(1); y=X(2); z=X(3); %

12、 % Y= X(4), X(7), X(10); % X(5), X(8), X(11); % X(6), X(9), X(12); % f=zeros(9,1); % f(1)=SIGMA*(y-x); f(2)=-x*z+R*x-y; f(3)=x*y-BETA*z; % % Jac=-SIGMA,SIGMA,0; R-z,-1,-x; y, x,-BETA; % % f(4:12)=Jac*Y; % % Run Lyapunov exponent calculation: % % T,Res=lyapunov(3,lorenz_ext,ode45,0,0.5,200,0 1 0,10);

13、 % % See files: lorenz_ext, run_lyap. % % - % Copyright (C) 2004, Govorukhin V.N. % This file is intended for use with MATLAB and was produced for MATDS-program % http:/www.math.rsu.ru/mexmat/kvm/matds/ % lyapunov.m is free software. lyapunov.m is distributed in the hope that it % will be useful, bu

14、t WITHOUT ANY WARRANTY. % % % n=number of nonlinear odes % n2=n*(n+1)=total number of odes % n1=n; n2=n1*(n1+1); % Number of steps nit = round(tend-tstart)/stept); % Memory allocation y=zeros(n2,1); cum=zeros(n1,1); y0=y; gsc=cum; znorm=cum; % Initial values y(1:n)=ystart(:); for i=1:n1 y(n1+1)*i)=1

15、.0; end; t=tstart; % Main loop TT,YY = feval(fcn_integrator,rhs_ext_fcn,t t+5000,y); y=YY(size(YY,1),:); fprintf(y value); for i=1:n1 fprintf( %10.6f,Y(i); end; fprintf(n); for ITERLYAP=1:nit % Solutuion of extended ODE system T,Y = feval(fcn_integrator,rhs_ext_fcn,t t+stept,y); %if (mod(ITERLYAP,io

16、utp)=0) % fprintf(y value); %for i=1:n1 % fprintf( %10.6f,Y(i); % end; % fprintf(n); %end; % for ii=1:500 % ddy=lorenz_ext(1,y); % Y=ddy*0.001+y; % y=Y; % end; t=t+stept; y=Y(size(Y,1),:); % y=Y; for i=1:n1 for j=1:n1 y0(n1*i+j)=y(n1*j+i); end; end; % % construct new orthonormal basis by gram-schmidt % znorm(1)=0.0; for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)2; end; znorm(1)=sqrt(znorm(1); for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end; for j=2:n1 for k