具有状态依赖时滞的 DDE
以下示例说明如何使用 ddesd
对具有状态依赖时滞的 DDE(时滞微分方程)方程组求解。Enright 和 Hayashi [1] 将此 DDE 方程组用作测试问题。
方程组为:
的历史解函数是解析解
方程中的时滞仅出现在 项中。时滞仅取决于第二个分量 的状态,因此这些方程构成状态依赖时滞方程组。
要在 MATLAB® 中求解此方程组,您需要先编写方程组、时滞和历史解的代码,然后再调用时滞微分方程求解器 ddesd
,该求解器适用于具有状态依赖时滞的方程组。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写时滞代码
首先,编写一个函数来定义方程组中的时滞。此方程组中唯一存在的时滞位于项 中。
function d = dely(t,y) d = exp(1 - y(2)); end
注意:所有函数都作为局部函数包含在示例的末尾。
编写方程代码
现在,创建一个函数来编写方程的代码。此函数应具有签名 dydt = ddefun(t,y,Z)
,其中:
t
是时间(自变量)。y
是解(因变量)。Z(n,j)
对时滞 求近似值,其中时滞 由dely(t,y)
的分量j
给出。
求解器会自动将这些输入传递给该函数,但是变量名称决定如何编写方程代码。在这种情况下:
Z(2,1)
function dydt = ddefun(t,y,Z) dydt = [y(2); -Z(2,1)*y(2)^2*exp(1 - y(2))]; end
编写历史解代码
接下来,创建一个函数来定义历史解。历史解是时间 的解。
function v = history(t) % history function for t < t0 v = [log(t); 1./t]; end
求解方程
最后,定义积分区间 并使用 ddesd
求解器对 DDE 求解。
tspan = [0.1 5]; sol = ddesd(@ddefun, @dely, @history, tspan);
对解进行绘图
解结构体 sol
具有字段 sol.x
和 sol.y
,这两个字段包含求解器在这些时间点所用的内部时间步和对应的解。(如果您需要在特定点的解,可以使用 deval
来计算在特定点的解。)
使用历史解函数绘制两个解分量对时间的图,以计算积分区间内的解析解来进行比较。
ta = linspace(0.1,5); ya = history(ta); plot(ta,ya,sol.x,sol.y,'o') legend('y_1 exact','y_2 exact','y_1 ddesd','y_2 ddesd') xlabel('Time t') ylabel('Solution y') title('D1 Problem of Enright and Hayashi')
局部函数
此处列出了 DDE 求解器 ddesd
为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。
function dydt = ddefun(t,y,Z) % equation being solved dydt = [y(2); -Z(2,1).*y(2)^2.*exp(1 - y(2))]; end %------------------------------------------- function d = dely(t,y) % delay for y d = exp(1 - y(2)); end %------------------------------------------- function v = history(t) % history function for t < t0 v = [log(t); 1./t]; end %-------------------------------------------
参考
[1] Enright, W.H. and H. Hayashi.“The Evaluation of Numerical Software for Delay Differential Equations.”In Proceedings of the IFIP TC2/WG2.5 working conference on Quality of numerical software: assessment and enhancement.(R.F. Boisvert, ed.).London, UK:Chapman & Hall, Ltd., pp. 179-193.
另请参阅
ddensd
| ddesd
| dde23
| deval