一阶时滞微分方程三种求解方法的MATLAB实现及稳定性分析

前言： 大学期间只学习过《常微分方程》，没想到有些学校竟然还学《时滞微分方程》，于是找到一本由内藤敏机(日本)等著，马万彪等译的《时滞微分方程——泛函数微分方程引论》(有需要的可以私聊，优快云貌似上传不了书籍，说侵权emmm)，看着头秃，不过受到不少启发，尤其是对Logistic方程的改进，真真是长见识了。没找到有人用欧拉法解一阶时滞微分方程的，于是一不做二不休便用MATLAB实现了一下下，顺便改进了点点。

一阶时滞微分方程在求解时常常会解着解着就变成了解超越方程了，而超越方程很难求解，很头疼，因此常常都是计算方程的数值解，而非解析解。本文也只计算数值解。

对于一阶时滞微分方程（方程(1)）
$$x^{\prime }=-px\left(t-\tau \right)+\frac{qx\left(t\right)}{r+x^{\alpha }\left(t\right)},t\ge 0,$$
这里$p,q,r$都是正的常数，$\tau$是时滞量，α∈N={1,2,3,...}\alpha \in N=\{1,2,3,...\}α∈N={1,2,3,...}，并且满足
$$\frac{q}{p}>r.$$
啊，不想打字了，直接放我写的报告的截图了。

取步长








当$p=10,q=2,r=0.1，\alpha =3,r=0.01$时，写出上述方程的欧拉法和改进欧拉法的Matlab程序，并进行稳定性分析和与解析解比较。

欧拉法是一种解决数值微分方程的最基本的方法，其基本思想是迭代。其中分为前进的欧拉法、后退的欧拉法、改进的欧拉法，本文使用的是前进的欧拉法。所谓迭代，就是逐次替代，最后求出所要求的解，并达到一定的精度。误差可以很容易地计算出来。

该方法有明显的优势：单步、显式、一阶求导精度、截断误差为二阶，以上优势能够大大减小运算的难度。但是不可避免的，欧拉法也存在着一定的缺陷。欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。

为了提高精度，需要在欧拉格式的基础上进行改进。本文采用区间两端的函数值的平均值作为直线方程的斜率得到改进的欧拉格式。改进欧拉法的精度为二阶。

除了改进的欧拉格式以外，还有许多比较好的改进方法，龙格库塔法便是其中一种。在数值分析中，龙格库塔法是用于模拟微分方程的解的重要的一类隐式或显式迭代法。实际上，龙格库塔法是欧拉法的一种推广，向前欧拉格式将导数项简单取为 ，而龙格库塔法则是在区间 内多取几点，将他们的斜率加权平均，作为导数的近似。

综合使用以上的欧拉法、改进欧拉法和龙格库塔法，取初值为1，步长 $h=0.5$，求解时滞微分方程(1)。编程时，我们先利用三中方法求解 $u_k$的值，由于$u(t)=x(\tau t)$ ，因此为了得到 $x(t)$，需要对解得的$u(t)$ 做如下逆变换
$$x(t)=u(\frac{t}{\tau })$$
得到的结果对比图如图1-4所示。

图1 方程(1)的解，步长$h=1/2$

图2 方程(1)的解，步长 $h=1/4$

图2 方程(1)的解，步长 $h=1/8$

图2 方程(1)的解，步长 $h=1/16$

对模型做稳定性分析，计算$\frac{q^2}{h\alpha \tau }$ 与$p(q-pr)-q^2$ 的值分别为266.67与6，因此该方程的解只有当 小于某个值时才会稳定，而$\tau =0.01$ 恰好在该范围内，因此解是稳定的。

①欧拉法代码：

function x = Eu(a,b,m,x0)
% a,b为x的取值范围
% m为单位时间的分段数，如将单位时间分为2段则m=2
% x0为x初始时刻的取值

%% 参数设置
p = 10;
q = 2;
r = 0.1;
alpha=3;
tao = 0.01;

h = 1/m;   % 步长
a_u = a/tao;   % 这是u的起始时间
b_u = b/tao;    % 这是u的终止时间
n = (b_u-a_u)/h;   % 总区间段数

t0 = a;   % 初始时间
u0 = x0;   % 初始值

%% 初始化0时刻之前的那些时间的值，有m个，再加上0时刻，总共有m+1个
u(1,1:m+1) = u0;     % 这是u_k
u(2,1:m+1) = -m:0;    % 这是t

%% 迭代计算后面时刻的值
for k = 0:n-1
    u(2,m+1+k+1) = k+1;   
    u_k = u(1,m+1+k);     % u_k
    u_km = u(1,m+1+k-m);   % u_{k-m}
    u(1,m+1+k+1) = (1 + q*h*tao/(r + (u_k)^alpha) )*u_k - tao*h*p*u_km;
end
u(2,:) = h*u(2,:);   % 转为真实时间

%% 转为x(t)
x(1,1) = x0;
x(2,1) = t0;
for k = 0:n-1
    x(1,k+2) = u(1,m+1+k+1);
    x(2,k+2) = u(2,m+1+k+1)*tao;
end
end

②改进欧拉法的代码：

function x = Eu_improve(a,b,m,x0)
% a,b为x的取值范围
% m为单位时间的分段数，如将单位时间分为2段则m=2
% x0为x初始时刻的取值

%% 参数设置
p = 10;
q = 2;
r = 0.1;
alpha=3;
tao = 0.01;

h = 1/m;    % 步长
a_u = a/tao;   % 这是u的起始时间
b_u = b/tao;    % 这是u的终止时间
n = (b_u-a_u)/h;   % 总区间段数

t0 = a;
u0 = x0;

%% 初始化0时刻之前的那些时间的值，有m个，再加上0时刻，总共有m+1个
u(1,1:m+1) = u0;     % 这是u_k
u(2,1:m+1) = -m:0;    % 这是t

%% 迭代计算后面时刻的值
for k = 0:n-1
    u(2,m+1+k+1) = k+1;   
    u_k = u(1,m+1+k);     % u_k
    u_km = u(1,m+1+k-m);   % u_{k-m}
    u_km1 = u(1,m+1+k-m+1);  % u_{k-m+1}
    a_k1 = q*tao*u_k/(r + (u_k)^alpha) - p*tao*u_km;
    a_k2 = q*tao*(u_k+h*a_k1)/( r + (u_k+h*a_k1)^alpha ) - p*tao*u_km1; 
    u(1,m+1+k+1) = u_k +h/2*(a_k1 + a_k2);
end
u(2,:) = h*u(2,:);   % 转为真实时间

%% 转为x(t)
x(1,1) = x0;
x(2,1) = t0;
for k = 0:n-1
    x(1,k+2) = u(1,m+1+k+1);
    x(2,k+2) = u(2,m+1+k+1)*tao;
end
end

③龙格库塔法

function x = LK(a,b,x0)

tao = 0.01;
lags = tao;
history = x0;
tspan = [a,b];
sol = dde23(@myddefun,lags,history,tspan);
x = [sol.y;sol.x];

% figure
% [t,x] = ode45(@myfun,[0 5],1);
% plot(t,x)
end


% function dxdt = myfun(t,x)
% p = 10;
% q = 2;
% r = 0.1;
% alpha=3;
% tao = 0.01;
% dxdt = -p*x.*(t-tao) + q*x./(r + (x).^alpha);
% 
% end

④主程序

clear,clc

% 参数设置
m =2;    % 步长h=1/m
x0 = 1;  % 初始值
a = 0;   % 区间左端点
b = 10;   % 区间右端点

% 欧拉法
x1 = Eu(a,b,m,x0);

% 改进欧拉法
x2 = Eu_improve(a,b,m,x0);

% 龙格-库塔法
x3 = LK(a,b,x0);

% 平稳性分析
p = 10;
q = 2;
r = 0.1;
alpha=3;
tao = 0.01;
h = 1/m;   % 步长
suan1 = q^2/(h*alpha*tao);
suan2 = p*(q-p*r)-q^2;
if suan1 < suan2
    disp('对任意tau>=0都是渐进稳定的')
else
    disp('在[0,tau_0)上是渐进稳定的')
end

% 绘图
plot(x1(2,:),x1(1,:),'*');
hold on
plot(x2(2,:),x2(1,:),'x');
plot(x3(2,:),x3(1,:),'o');
legend('欧拉法','改进欧拉法','龙格-库塔法','location','best')
xlabel('t')
ylabel('x(t)')
set(gca,'linewidth',1.1)

(4)比较两种方法的差别。
通过综合观察图1-4，不难发现，对于一阶时滞微分方程(1)的求解，两种方法得到的结果几乎一致，且都与准确解吻合得非常好。因此可以认为，即使在大多数情况下或者说从理论的角度来讲，欧拉法存在着明显的不足，但是在某些情况下欧拉法求得的解仍然比较好，且由于欧拉法相较于改进欧拉法以及龙格库塔法而言逻辑简单通俗易懂且容易实现的特点，往往更具有实用价值。