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)


相关推荐

  • c语言之异或运算_c语言运算符优先级表

    c语言之异或运算_c语言运算符优先级表c语言之异或运算异或运算,计算机相关专业比较熟悉了。相同为0,不同为1.结合计算机内部的位运算,a^a=0;与本身异或是为0的。有关的知识运用到数据交互中去。voidint_swap(int*x,int*y){ *y=*x^*y;//step1 *x=*x^*y;//step2 *y=*x^*y;//step3}运用这个函数就能完成两个数据交换。但是并没有提高时间复杂度和空间复杂度,有关书籍上称之为智力游戏。我们来看看数据的变化。假设*x=a,*y=b.*x*y

    2022年10月31日
  • 负载均衡、Apache 负载均衡配置[通俗易懂]

    负载均衡、Apache 负载均衡配置[通俗易懂][1]Apache负载均衡设置方法mod_proxy使用介绍一般来说,负载均衡就是将客户端的请求分流给后端的各个真实服务器,达到负载均衡的目的。还有一种方式是用两台服务器,一台作为主服务器(Master),另一台作为热备份(HotStandby),请求全部分给主服务器,在主服务器当机时,立即切换到备份服务器,以提高系统的整体可 第一次看到这个标题时我也很惊讶,Apache居然还能做负载

  • lena图像下载「建议收藏」

    lena图像下载「建议收藏」 http://www.ece.rice.edu/~wakin/images/

  • 数据接口-免费版(股票数据API)「建议收藏」

    获取股票数据的源头主要有:数据超市、雅虎、新浪、Google、和讯、搜狐、ChinaStockWebService、东方财富客户端、证券之星、网易财经。数据超市2016年5月6日更新。根据最近频繁出现的数据超市,可以无限制获取相关数据,而不再需要使用爬虫等方式获取,这样不仅节省了极大资源,也有利于遍历数据。具体的方法不再赘述,列出来相关网站清单,开发者可自行到这些网站查询调用方法。聚合数据 htt…

  • 程序员去外包公司有前途吗_程序员去外包是不是就废了

    程序员去外包公司有前途吗_程序员去外包是不是就废了虽然大部分人都抵制外包,但是很多人,尤其是萌新,并不清楚外包的主要缺点。我这里简单说一下。程序员去外包公司有前途吗?不能说去了外包公司就完全没有前途了,主要看个人能力,外包的工作内容,大多十分碎片化,甚至是机械化。因为如果这个工作内容真的很完整、成块儿,那正式工就做掉了。正式工做掉的理由有两个:完整工作内容有利于他,去构建业务认知。完整内容拆分出来外包,需要进行进行大量的沟通与团队协作,不利于整体效率。那么有没有办法避免碎片化呢?答案是有的。一方面可以表现出自身能力,获取正式团队

  • 传感器尺寸

    传感器尺寸传感器尺寸   图像传感器感光区域的面积大小。这个尺寸直接决定了   整个系统的物理放大率。如:1/3”、1/2”等。绝大多数   模拟相机的传感器的长宽比例是4:3(H:V),数字相机   的长宽比例则包括多种:1:1,16:9,3:2etc。   注意相机尺寸的1”=16mm≠25.4mm…

发表回复

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

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