matlab求解高阶常微分方程_状态依赖时滞微分方程的动力学研究

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

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE稳定放心使用

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

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

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

取步长
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
p = 10 , q = 2 , r = 0.1 , α = 3 , r = 0.01 p=10,q=2,r=0.1,\alpha=3,r=0.01 p=10,q=2,r=0.1α=3,r=0.01时,写出上述方程的欧拉法和改进欧拉法的Matlab程序,并进行稳定性分析和与解析解比较。

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

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

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

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

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

在这里插入图片描述
图2 方程(1)的解,步长 h = 1 / 16 h=1/16 h=1/16

对模型做稳定性分析,计算 q 2 h α τ \frac{q^2}{h\alpha \tau} hατq2 p ( q − p r ) − q 2 p(q-pr)-q^2 p(qpr)q2 的值分别为266.67与6,因此该方程的解只有当 小于某个值时才会稳定,而 τ = 0.01 \tau=0.01 τ=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)的求解,两种方法得到的结果几乎一致,且都与准确解吻合得非常好。因此可以认为,即使在大多数情况下或者说从理论的角度来讲,欧拉法存在着明显的不足,但是在某些情况下欧拉法求得的解仍然比较好,且由于欧拉法相较于改进欧拉法以及龙格库塔法而言逻辑简单通俗易懂且容易实现的特点,往往更具有实用价值。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/187202.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • A Painless Q-learning Tutorial (一个 Q-learning 算法的简明教程)

    本文是对 http://mnemstudio.org/path-finding-q-learning-tutorial.htm 的翻译,共分两部分,第一部分为中文翻译,第二部分为英文原文。翻译时为方便读者理解,有些地方采用了意译的方式,此外,原文中有几处笔误,在翻译时已进行了更正。这篇教程通俗易懂,是学习理解Q-learning算法工作原理的绝佳入门材料。

  • 科学组合,系统学习

    科学组合,系统学习

  • 图片管理系统–CDN源站图片管理

    这是Java开发的跨域的图片上传、图片删除解决办法;并提供了5个针对商品的图片管理功能:浏览、检查、切割、删除、上传。对于图片上传、图片删除是没有权限控制的;对于5个图片管理功能演示,使用JasigCAS进行单点登录。技术框架:SpringMVC下载连接一:图片管理系统源码下载地址:http://pan.baidu.com/s/11

  • 小波变换对图像的分解与重构(含matlab代码)

    小波变换对图像的分解与重构(含matlab代码)01小波变换原理所谓的小波的小是针对傅里叶波而言,傅里叶波指的是在时域空间无穷震荡的正弦(或余弦波)。相对而言,小波指的是一种能量在时域非常集中的波,它的能量有限,都集中在某一点附近,而且积分的值为零,这说明它与傅里叶波一样是正交波。举…

  • charles导致mac无法上网_使用不同的MAC地址上网

    charles导致mac无法上网_使用不同的MAC地址上网前言charles关闭后,发现网页突然打开了,那大概率是设置了代理,但明明已经关闭了charles,这是由于mac网络偏好设置中,使用的是手动代理,将其改为自动即可解决方法1打开网络偏好设置,

  • 日期函数months_between的用法[通俗易懂]

    日期函数months_between的用法[通俗易懂]MONTHS_BETWEEN(date1,date2)用于计算date1和date2之间有几个月。如果date1在日历中比date2晚,那么MONTHS_BETWEEN()就返回一个正数。如果date1在日历中比date2早,那么MONTHS_BETWEEN()就返回一个负数。如果date1和date2日期一样,那么MONTHS_BET…

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号