MatLab求解延迟微分方程的注意事项.doc

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

2、有关的延迟微分方程；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就会找看到函数作者

4、所写的一个入门例子：function 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)

5、 = y_1(t-1)+y_2(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 tau=1;%给定延迟时间 history=0;%初始值 tapan0,10;%求解时间范围 sol=dde23(yanchi,tau,history,tapan);%延迟问题求解 plot(sol.x,sol.y)；%作图下面附上了图片 x(0)=0 和 x(0)=2 的情况显然初始值不同结果不同，这就是为什么需要初始值的情况延迟微分方程是：dx(t)/dt=0.2

6、*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 - Lyapunov 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_ex

7、t_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 Jacobian 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

8、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); % % 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); % % See files: lorenz_ext, run_lyap.