延迟函数

使用MATLAB求解延时微分方程的两种：DDE23和SimuLink有些不同点需要注意，否则结果会出现错误   使用MATLAB来求解延迟微分方程是在生物数学和化学计算求解中经常遇到的事，在其它领域也比较常见。我所知道的，在MATLAB中求解微分方程有三种方法：1.使用ode45（龙格-库塔法的一个变种）求解，这时用一个数组，y的延迟项，但是c的值要考虑步长，再代入方程就能实现延时效应；2.使用dde23求解常数延时方程、使用ddesd可以用来求解延迟与时间t有关的延迟微分方程；3.使用SimuLink建模求解，SimuLink是求解广义微分代数系统的通用工具，功能很强大，但是看惯了编程指令的人可能不大习惯，调试似乎也不太方便。        既然本文专门讨论求解延迟微分方程，就先介绍一下专用工具dde23，dde系列求解函数是由Southern Methodist University 的L.F. Shampine 和S. Thompson根据他们早期使用Fortran编写的Fortran 90 DDE Solver 移植到MATLAB上的，从MATLAB6.5开始加入MATLAB的官方发行版，根据薛定宇教授在其几本关于MATLAB的著作中提到的，该函数返回的sol中结构体sol.x和sol.y均为按行排列，与ode45等不同，不太规范（没办法，因为这个函数本来就不是Mathworks的官方作品），不过这一点已经不大可能得到改进了，因为L.F. Shampine 和S. Thompson已经决定停止维护这个文件。如果您想进一步了解该函数，可以访问它的主页。     MATLAB帮助中关于该函数的介绍不很清楚，如果需要进一步了解这个函数，需要下载作者为其写的手册。下面以MATLAB中所附的一个例程来说明这个函数与Simulink建模求解的不同。    在MATLAB prompt中键入edit ddex1就会找看到函数作者所写的一个入门例子：   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) = 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 <= 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 is 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],@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); 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<0任意时刻的值。也就是说，仿真开始时刻的 这样，仿真从0时刻开始y_1就以初始为1的变化率上升。同理，y_2的情况也一样。        而在simulink中，我们搭建等效模型如图：   建模方程是正确的，各积分模块的初始值为1，延迟模块的延时值设为1和0.2，但是得出的结果与DDE23有非常大的差别。经过数次更正求解器参数依然没有改观。经过对比仿真结果可以发现：不同于DDE23，SimuLink的仿真时间设为T=0，此时，y_1的初始值为1，但是Y_history的值没有得到体现，也就是说y'_1的值为0，只有到T=1时刻以后，才会有y'_1(t)=y(t-1)，方程延时效应才正式得到体现。这与DDE23的基本设定是有差别的。 如果方程模型中只有一个延时常数1，将simulink的仿真时间定为0-6秒，其中1-6秒与DDE23的仿真结果相同。在本例中有两个时间常数，这样做行不通。那么我们可以通过设定两个延时模块的initial output选项分别为y_1(-1)=1，y_2(-0.2)=1，这样就会得到与DDE23相同的结果。   因此，在使用MATLAB求解延迟微分方程时，应该先弄清楚延迟方程的具体意义，再选择适当的求解方式，才能得到正确的结果。 顺便提一下，MATLAB被广为诟病的除了昂贵的价格之外，最突出的就是透明度较低的内部算法，使用者往往搞不清算法的内部结构（很多算法甚至不予公开），也就无从判断求解结果的正确性，这时，转向一些开源软件（SciLab、MAXIMA、Octave、SAGE）和开放式求解库可以解决问。