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

相关推荐

  • 彻底弄懂C语言数组名

    彻底弄懂C语言数组名先定义一个一维数组inta[]={0,1,2,3,4,5,6,7,8,9};一、数组名是什么数组名的值是数组首元素的指针常量。数组名不是指针,但大多数使用到数组名的地方,编译器都会把数组名隐式转换成一个指向数组首元素的指针来处理。只有两种情况下例外:第一种是对数组名使用sizeof运算符sizeof(a)这将会得到整个数组所占的内存大小,a是长度为10的int(4字节

  • touchstart与click同时触发

    touchstart与click同时触发产生冲突的原因我们可以给某个元素同时绑定touchstart和click事件,但这将会导致本篇文章解决的问题–这两个事件在移动设备上会发生冲突。由于移动设备能够同时识别touchstart和click事件,因此当用户点击目标元素时,绑定在目标元素上的touchstart事件与click事件(约300ms后)会依次被触发,也就是说,我们所绑定的回调函数会被执行两次!。…

  • html实现滑动解锁_js滑动解锁(原创)

    html实现滑动解锁_js滑动解锁(原创)varbox=document.querySelector(“#box”);varcontent=document.querySelector(“#content”);varshadow=document.querySelector(“#shadow”);vartip=document.querySelector(“#tip”);varblock=document.q…

  • asp.net中DropDownList控件各种属性研究汇总

    asp.net中DropDownList控件各种属性研究汇总.aspx代码如下:AutoPostBack=”True”>AutoPostBack=”true”onselectedindexchanged=”DropDownList2_SelectedIndexChanged”>

    2022年10月17日
  • Vim如何撤销和回退[通俗易懂]

    Vim如何撤销和回退[通俗易懂]vim中的撤销与回退u,撤销ctrl+R,回退(rollback,回滚、回退)

  • Java分布式锁

    Java分布式锁Java分布式锁我的理解应该叫集群锁或者跨实例锁锁的作用是在多线程情况下,控制线程同步访问变量,执行代码块、方法,例如synchronized,在单个jvm进程中,这样是奏效的。但是在分布式环境中,单个服务往往都是要部署多台实例的,在有多个jvm进程的集群里,synchronized就达不到我们的要求了。synchronized只能控制当前jvm进程中的线程,对于其它jvm进程中的线程,它无能为力。也就是说有可能一个jvm中的线程是同步执行的,在此过程中,或许会有集群里其它jvm的线程执行到

发表回复

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

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