营销型网站建设要点,网站整合营销推广,智能建站开发,设计免费素材网站有哪些2023年数学建模美赛A题#xff08;A drought stricken plant communities#xff09;分析与编程 2023年数学建模美赛D题#xff08;Prioritizing the UN Sustainability Goals#xff09;分析与编程 特别提示#xff1a; 1 本文介绍2023年美赛题目#xff0c;进行深入分析… 2023年数学建模美赛A题A drought stricken plant communities分析与编程 2023年数学建模美赛D题Prioritizing the UN Sustainability Goals分析与编程 特别提示 1 本文介绍2023年美赛题目进行深入分析 2 本文首先对 A题进行深入分析其它题目分析详见专题讨论 3 最新更新增加了多种群LK模型和例程。 文章目录2023年数学建模美赛A题A drought stricken plant communities分析与编程1. A题A drought stricken plant communities遭受旱灾的植物群落1.1 背景1.2 要求2. 问题分析3. 种间竞争关系理论 Lotka-Volterra模型3.1 Malthus人口理论3.2 捕食者-猎物系统3.3 Lotka-Volterra模型3.4 Lotka-Volterra模型的数值模拟3.4 三种群 Lotka-Volterra模型的数值模拟4. Lotka-Volterra模型与全球植被模式4.1 现存生物量一般动力学模型4.2 包含n个物种的全球植被模式的一般动力学方程5. 参考数据6. 参考资料3. 微分方程模型求解3.1 基本概念3.2 微分方程的数学建模3.3 微分方程的数值解法4. SciPy 求解常微分方程组4.1 一阶常微分方程组模型4.2 scipy.integrate.odeint() 函数5. 求洛伦兹Lorenz方程的数值解6.1 例题 2求洛伦兹Lorenz方程的数值解6.2 洛伦兹Lorenz方程问题的编程步骤6.3 洛伦兹Lorenz方程问题 Python 例程6.4 洛伦兹Lorenz方程问题 Python 例程运行结果2023年数学建模美赛A题A drought stricken plant communities分析与编程 1. A题A drought stricken plant communities遭受旱灾的植物群落
1.1 背景
不同的植物群落对压力的反应不同。例如草原对干旱非常敏感。干旱发生的频率和严重程度不同。但大量的观察表明不同物种的数量对植物群落如何在连续多代干旱周期中的适应能力起到了重要作用。在一些仅有单一物种的植物群落中后代不像有 4种或更多物种的群落中的个体植物那样容易适应干旱条件。这些观察引出了许多问题例如对于一个植物群落要从这种局部生物多样性中受益最少需要多少种物种随着物种数量的增加这种现象如何发展这对植物群落的长期生存性意味着什么 1.2 要求
考虑到植物群落中干旱适应性与物种数量的关系您的任务是探索和更好地理解这一现象。具体地您应该 开发一个数学模型预测植物群落随着不同的不规则天气周期的变化。在降水充足的时期应包括降雨的时间。该模型应考虑干旱周期中不同物种之间的相互作用。 探讨你能从你的模型中得出什么结论关于植物群体与更大环境的长期相互作用。考虑以下问题 要使植物群落受益需要的不同植物物种数量是多少随着物种数量的增加会发生什么社区中的物种类型如何影响你的结果未来天气周期中干旱发生的频率和变化范围的影响是什么如果干旱较少物种数量对总人口的影响是否相同污染和栖息地减少等其他因素如何影响你的结论您的模型表明应该采取什么措施以确保植物群落的长期生存力对更大环境的影响是什么 2. 问题分析 这是一道微分方程建模题目建立模型是关键模型求解并不难基于模型的分析和讨论可以发挥想象力。需要先找到相关研究论文根据论文中提出物种与环境的关系的原理模型建立微分方程的数学模型。论文中会给出具体的数学模型可能是偏微分方程能够求解就直接用如果不会就简化为常微分方程也可以。微分方程是描述系统的状态随时间和空间演化的数学工具。本题显然是研究几种物种的数量随时间的变化规律。**特别注意**给定初始条件的一阶常微分方程组的标准形式是 {dydtf(y,t)y(t0)y0\begin{cases} \begin{aligned} \frac{dy}{dt} f(y,t)\\ y(t_0) y_0 \end{aligned} \end{cases} ⎩⎨⎧dtdyf(y,t)y(t0)y0 微分方程是微分方程组式中的 y 是数组向量有几个物种就有几个变量yi(t)y_i(t)yi(t) 表示物种 i 的总量随时间 t 的变化。 也就是说可以先建立一种或两种物种的模型分析变化趋势再依次增多物种数量分析变化趋势。**特别注意**构造外部条件及降雨量随时间的变化函数。 降雨量函数要满足题目要求1不同的不规则天气周期既要有干旱周期又要有降水周期还要有间隔周期2干旱发生的频率和变化范围。 如何构造合理的降雨量函数可以体现能力和创新的。思路一是设计的降雨量函数包含不同的可能模式思路二是降雨量函数包括地球上典型的干旱模式思路三是找到非洲地区的降雨量统计。**特别注意**污染和栖息地减少等其他因素的影响需要在模型中增加一个系数或者一项反映污染和栖息地减少的影响。直接在模型中增加一个系数比较简单。 进一步地可以构思污染与种群总数相关是一个随种群总数变化的系数。**特别注意**采取什么措施以确保植物群落的长期生存力应该通过模型研究得到结论。 建议首先查找资料找到一种或多种靠谱的措施然后将其量化为一个系数或一项加入模型跟污染系数的原理是一致的通过模型研究措施的影响证明采取的措施是有效的。 简单地既然污染会破坏环境那么防止污染就可以保护环境可以抑制污染系数的增大。 又如水土保持能否量化为模型参数食物链能否构造模型微分方程模型求解详见本文后续章节及博客 Python小白的数学建模课-09.微分方程模型https://youcans.blog.csdn.net/article/details/117702996Python小白的数学建模课-11.偏微分方程数值解法https://youcans.blog.csdn.net/article/details/119755450 本题要研究时间的变化因此跟微分方程边值问题没什么关系。 3. 种间竞争关系理论 Lotka-Volterra模型
常微分方程与生态学有什么关联生态学理论一直以种群动态为根基。它涌现出所谓生物多样性、空间分布格局、种间互作也是功能性状、物候等话题在应用生态学中的归宿。
3.1 Malthus人口理论
1798年Malthus提出的人口理论被认为是生态学第一定律。Malthus指出一切生物在“不受外力”的影响下都以恒定的速率 k 增长。用微分方程可表达为 dxdtkx,k0\frac{dx}{dt} kx, k0 dtdxkx,k0
该方程可直接用定积分或不定积分求解结果为 xek(t−t0)x0xe^{k(t-t_0)}x_0 xek(t−t0)x0 这是一个指数曲线 称为**「自然繁殖normal reproduction** 。
由于资源有限当种群个体数过大时竞争激烈故种群增长率会有所下降。所以 k 不是一个恒值而是关于 x 的函数 kf(x)kf(x) kf(x)
函数f(x)应当是在x0范围内随x增加而单调递减的函数。根据麦克劳林展开式
f(x)f(0)df(x)dx∣x0(x)o(x)≈a−bxf(x)f(0)\frac{df(x)}{dx}|_{x0} (x)o(x)\approx a-bx f(x)f(0)dxdf(x)∣x0(x)o(x)≈a−bx
当x足够小时任何光滑函数可用线性函数作近似因此:
dxdt(a−bx)x\frac{dx}{dt}(a-bx)x dtdx(a−bx)x
求解微分方程可得
xae−atbx\frac{a}{e^{-at}b} xe−atba 3.2 捕食者-猎物系统
Lotka1910和Volterra1926分别提出了捕食者–猎物系统的种群动态模型称为Lotka–Volterra模型。该模型假设捕食者、猎物相遇的几率与双方的种群个体数成正比捕食者捕食猎物的频次与相遇频次成正比捕食者种群增长率与捕食猎物的频次成正比捕食者以恒定的速率死亡。 dxdtkx−axydydt−lybxy\frac{dx}{dt}kx-axy\\ \frac{dy}{dt}-lybxy dtdxkx−axydtdy−lybxy
这个二元系统捕食者猎物动态会如何变化
► 捕食者、猎物种群动态的相位曲线如何
► 捕食者、猎物种群动态是达到稳定平衡还是有限环还是混沌或者说相位曲线会呈螺旋形造成系统崩溃或收敛到一个稳定点
► 受到外界扰动时捕食者、猎物种群动态还能保持稳定吗
1 为什么到达平衡的时间经常是无穷大
2 初值对Lotka–Volterra模型行为的影响
3 无扰动下Lotka–Volterra有限环大小不变
3.3 Lotka-Volterra模型
Lotka-Volterra模型Lotka-Volterra种间竞争模型是logistic模型阻滞增长模型的延伸。现设定如下参数
N1、N2分别为两个物种的种群数量K1、K2分别为两个物种的环境容纳量r1、r2 分别为两个物种的种群增长率
依逻辑斯蒂模型有如下关系 dN1dtr1∗N1∗(1−N1K1)\frac{dN_1}{dt}r_1 * N_1 * (1 - \frac{N_1}{K_1}) dtdN1r1∗N1∗(1−K1N1)
其中N/K 可以理解为已经利用的空间称为“已利用空间项”则 (1−N/K)(1-N/K)(1−N/K) 可以理解为尚未利用的空间称为“未利用空间项”。 当两个物种竞争或者利用同一空间时“已利用空间项”还应该加上N2种群对空间的占用。则 KaTeX parse error: Undefined control sequence: \alfa at position 53: …ac{N_1}{K_1} - \̲a̲l̲f̲a̲ ̲\frac{N_2}{K_2}…
其中α物种2对物种1的竞争系数即每个N2个体所占用的空间相当于α个N1个体所占用空间。 则有β物种1对物种2的竞争系数即每个N1个体所占用的空间相当于β个N2个体所占用空间。则另有
dN2dtr2∗N2∗(1−N2K2−βN1K1)\frac{dN_2}{dt}r_2 * N_2 * (1 - \frac{N_2}{K_2} - \beta \frac{N_1}{K_1} ) dtdN2r2∗N2∗(1−K2N2−βK1N1)
当物种N1种群物种1的环境容纳量为K1时N1种群中每个个体对自身种群的增长抑制作用为1/K1 同理N2种群中每个个体对自身种群的增长抑制作用为1/K2。 另外从1、2两个方程以及α、β的定义中可知
N2种群中每个个体对N1种群的影响为α/K1N1种群中每个个体对N2种群的影响为β/K2 因此当物种2可以抑制物种1时可以认为物种2对物种1的影响 物种2对自身的影响即 α/K1 1/K2。 整理后得K2 K1/α同理有物种2不能抑制物种1K2 K1/α物种1可以抑制物种2K1 K2/β物种1不能抑制物种2K1 K2/β
在竞争的过程中由于K1、K2、α 以及 β 的数值不同可能会产生如下四种结果 1物种1能抑制物种2物种2能抑制物种1两物种都有可能得胜 2物种1不能抑制物种2物种2能抑制物种1物种 2 得胜 3物种1能抑制物种2物种2不能抑制物种1物种 2 得胜 4物种1不能抑制物种2物种2不能抑制物种1两物种都不能抑制对方形成稳定平衡。 结果4是一个稳定的平衡无论N1和N2种群数量的组合(N1,N2)落在直角坐标系内哪一区域最终都将使得N1种群和N2种群的数量趋向平衡点。 3.4 Lotka-Volterra模型的数值模拟
引入如下的相对种群规模变量
$$ u\frac{U}{U^∗}\frac{eγ}{β}U\
v\frac{V}{V^∗}\frac{γ}{a}V\ $$
以及无量纲时间 ταβtταβtταβt
原来的 Lotka-Volterra 方程组可以被改写成
{dudτr∗(u−uv)dvdτ1r(uv−v)\begin{cases} \frac{du}{dτ} r*(u−uv)\\ \frac{dv}{dτ}\frac{1}{r}(uv−v) \end{cases} {dτdur∗(u−uv)dτdvr1(uv−v)
这里只剩下了一个自由的无量纲参数 rαβ。如此便说明在数值模拟时我们只需要调节1个而不是4个参数问题可以大大简化。r 值的大小和系统初始条件就可以完全确定解的行为。
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np # 导入 numpy包
import matplotlib.pyplot as plt # 导入 matplotlib包def dyLV(y, t, r,mu): # SIR 模型导数函数u, v ydu_dt r * (u- u*v) # du/dt r*(u-uv)dv_dt (u*v-v) / r # dv/dt (uv-v)/rreturn np.array([du_dt,dv_dt])# 设置模型参数
r [0.25, 0.5, 1.0, 2.0, 4.0]
mu 1.0
colorlist [green,cyan,blue,orange,red]
tEnd 20 # 预测日期长度
t np.arange(0.0, tEnd, 1) # (start,stop,step)
u0 1.0 # 初值
v0 1.2 # 初值
Y0 (u0, v0) # 微分方程组的初值plt.figure(figsize(10, 4))
for k in range(len(r)):yt odeint(dyLV, Y0, t, args(r[k],mu)) # SIS 模型plt.plot(t, yt[:, 0], colorcolorlist[k], labelrr {}.format(r[k]))plt.plot(t, yt[:, 1], --, colorcolorlist[k], labelr {}.format(r[k]))plt.xlabel(t)
plt.legend(locbest) # youcans
plt.show()下面几张图给出了不同 r 的取值对应的猎物和猎食者数量随时间的变化 3.4 三种群 Lotka-Volterra模型的数值模拟
{du1dtu1∗(1−u1−αu2−βu3)du2dtu2∗(1−βu1−u2−αu3)du3dtu3∗(1−αu1−βu2−u3)\begin{cases}\frac{du_1}{dt} u_1 * (1 - u_1 -\alpha u_2 -\beta u_3)\\\frac{du_2}{dt} u_2 * (1 - \beta u_1 - u_2 -\alpha u_3)\\\frac{du_3}{dt} u_3 * (1 - \alpha u_1 - \beta u_2 - u_3)\end{cases}\\ ⎩⎨⎧dtdu1u1∗(1−u1−αu2−βu3)dtdu2u2∗(1−βu1−u2−αu3)dtdu3u3∗(1−αu1−βu2−u3)
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np # 导入 numpy包
import matplotlib.pyplot as plt # 导入 matplotlib包def dyLV3(y, t, alfa, beta): # 3物种LK模型导数函数u1, u2, u3 ydu1_dt u1 * (1 - u1 - alfa*u2 - beta*u3)du2_dt u2 * (1- beta*u1 -u2 - alfa*u3)du3_dt u3 * (1- alfa*u1 -beta*u2 - u3)return np.array([du1_dt, du2_dt, du3_dt])# 设置模型参数
alfa, beta 1.5, -1.0
tEnd 10 # 预测长度
t np.arange(0.0, tEnd, 0.1) # (start,stop,step)
u10, u20, u30 0.2, 0.3, 0.5 # 初值
Y0 (u10, u20, u30) # 微分方程组的初值plt.figure(figsize(9,6))
# plt.subplot(121), plt.title(1. u(r))
yt odeint(dyLV3, Y0, t, args(alfa, beta)) # SIS 模型
plt.plot(t, yt[:,0], labelu1(t))
plt.plot(t, yt[:,1], labelu2(t))
plt.plot(t, yt[:,2], labelu3(t))
plt.xlabel(t)
plt.legend(locbest)
plt.show()4. Lotka-Volterra模型与全球植被模式
Svirezhev YLotka-Volterra模型与全球植被模式。Ecological Modelling2000135
4.1 现存生物量一般动力学模型
《资源环境数学模型手册》 (2019). 标准Lotka-Volterra模型
dBdtP−mR∗B\frac{dB}{dt} P - m_R*B dtdBP−mR∗B
式中B 为现存生物量P 是净生产率rR1/mRr_R1/m_RrR1/mR 是碳在生物量中实际停留时间一般来说PP(T,W,I,C)PP(T,W,I,C)PP(T,W,I,C)其中T是时间W是土壤含水量I是光照C是大气中碳浓度。
mRmG(B,Bs)m_R m G(B,B_s) mRmG(B,Bs)
式中BsB_sBs 为其他竞争物种的生物量可理解为某一物种的平均生理寿命。
4.2 包含n个物种的全球植被模式的一般动力学方程
dBidtPi(Bi,T,H)−miR∗Bi\frac{dB_i}{dt} P_i(B_i,T,H) - m_i^R*B_i dtdBiPi(Bi,T,H)−miR∗Bi
式中T 和H是这一空间位置的年平均温度和年降水量P_i 净初级生产力miR为第i类物种的平均生理寿命。
假设竞争系数是空间无关的即便可得到标准的Lotka-Volterra模型
∂Bi∂tBi∗[gi(T,H)−mi−∑γijBj\frac{\partial B_i}{\partial t} B_i * [g_i(T,H) - m_i - \sum {\gamma _{ij} B_j} ∂t∂BiBi∗[gi(T,H)−mi−∑γijBj 5. 参考数据 6. 参考资料
(1) 植被变化与降水量、降水变率的关系
(2) 长期降水量变化下荒漠草原植物生物量、多样性及其影响因素研究
(3) 荒漠植物幼苗对模拟降水量变化的响应
(4) 中亚地区1982年至2002年植被指数与气温和降水的相关性分析
(5) 内蒙古典型草原生长季内植物生长动态的数学模型与计算机模拟研究
(6) 降水量的季节分配对羊草草原群落地上部生物量影响的数学模型
(7) 干旱区荒漠植被生态需水量计算方法研究
(8) 干旱半干旱区植被生态需水量计算方法评述
(9) 基于数学模型的干旱区植被耗水量估算
(10) 乌审旗植被覆盖度动态变化及其与降水量的关系 (11) 归一化植被指数与降水量,土壤湿度的关系 以下内容只供参考
3. 微分方程模型求解
3.1 基本概念
微分方程是描述系统的状态随时间和空间演化的数学工具。物理中许多涉及变力的运动学、动力学问题如空气的阻力为速度函数的落体运动等问题很多可以用微分方程求解。微分方程在化学、工程学、经济学和人口统计等领域也有广泛应用。
具体来说微分方程是指含有未知函数及其导数的关系式。
微分方程按自变量个数分为只有一个自变量的常微分方程Ordinary Differential Equations和包含两个或两个以上独立变量的偏微分方程Partial Differential Equations。微分方程按阶数分为一阶、二阶、高阶微分方程的阶数取决于方程中最高次导数的阶数。微分方程还可以分为非齐次常变系数非线性初值问题/边界问题… 3.2 微分方程的数学建模
微分方程的数学建模其实并不复杂基本过程就是分析题目属于哪一类问题、可以选择什么微分方程模型然后如何使用现有的微分方程模型建模。
在数学、力学、物理、化学等各个学科领域的课程中针对该学科的各种问题都会建立适当的数学模型。在中学课程中各学科的数学模型主要是线性或非线性方程而在大学物理和各专业的课程中越来越多地出现用微分方程描述的数学模型。
数学建模中的微分方程问题通常还是这些专业课程中相对简单的模型专业课程的教材在介绍一个模型时往往都做了非常详细的讲解。只要搞清楚问题的类型、选择好数学模型建模和求解并不是很难而且在撰写论文时对问题背景、使用范围、假设条件、求解过程有大量现成的内容可以复制参考。
希望你在学习本系列之后会发现微分方程模型是数学建模中最容易的题型模型找教材建模找例题求解有例程讨论有套路论文够档次。 3.3 微分方程的数值解法
在学习专业课程时经常会推导和求解微分方程的解析解小白对微分方程模型的恐惧就是从高等数学“微分方程”开始经过专业课的不断强化而形成的。实际上只有很少的微分方程可以解析求解大多数的微分方程只能采用数值方法进行求解。
微分方程的数值求解是先把时间和空间离散化然后将微分化为差分建立递推关系然后反复进行迭代计算得到任意时间和空间的值。
如果你还是觉得头晕目眩我们可以说的更简单一些。建模就是把专业课教材上的公式抄下来求解就是把公式的参数输入到 Python 函数中。
我们先说求解。求解常微分方程的基本方法有欧拉法、龙格库塔法等可以详见各种教材撰写数模竞赛论文时还是可以抄几段的。本文沿用“编程方案”的概念不涉及这些算法的具体内容只探讨如何使用 Python 的工具包、库函数零基础求解微分方程模型。
我们的选择是 Python 常用工具包三剑客Scipy、Numpy 和 Matplotlib
Scipy 是 Python 算法库和数学工具包包括最优化、线性代数、积分、插值、特殊函数、傅里叶变换、信号和图像处理、常微分方程求解等模块。有人介绍 Scipy 就是 Python 语言的 Matlab所以大部分数学建模问题都可以用它搞定。Numpy 提供了高维数组的实现与计算的功能如线性代数运算、傅里叶变换及随机数生成另外还提供了与 C/C 等语言的集成工具。Matplotlib 是可视化工具包可以方便地绘制各种数据可视化图表如折线图、散点图、直方图、条形图、箱形图、饼图、三维图等等。
顺便说一句还有一个 Python 符号运算工具包 SymPy以解析方式求解积分、微分方程也就是说给出的结果是微分方程的解析解表达式。很牛但只能求解有解析解的微分方程所以你知道就可以了。 4. SciPy 求解常微分方程组
4.1 一阶常微分方程组模型
给定初始条件的一阶常微分方程组的标准形式是 {dydtf(y,t)y(t0)y0\begin{cases} \begin{aligned} \frac{dy}{dt} f(y,t)\\ y(t_0) y_0 \end{aligned} \end{cases} ⎩⎨⎧dtdyf(y,t)y(t0)y0
式中的 y 在常微分方程中是标量在常微分方程组中是数组向量。 4.2 scipy.integrate.odeint() 函数
SciPy 提供了两种方式求解常微分方程基于 odeint 函数的 API 比较简单易学基于 ode 类的面向对象的 API 更加灵活。
**scipy.integrate.odeint() **是求解微分方程的具体方法通过数值积分来求解常微分方程组。在 odeint 函数内部使用 FORTRAN 库 odepack 中的 lsoda可以求解一阶刚性系统和非刚性系统的初值问题。官网介绍详见 scipy.integrate.odeint — SciPy v1.6.3 Reference Guide 。
scipy.integrate.odeint(func, y0, t, args(), DfunNone, col_deriv0, full_output0, mlNone, muNone, rtolNone, atolNone, tcritNone, h00.0, hmax0.0, hmin0.0, ixpr0, mxstep0, mxhnil0, mxordn12, mxords5, printmessg0, tfirstFalse)odeint 的主要参数
求解标准形式的微分方程组主要使用前三个参数
func: callable(y, t, …) 导数函数 f(y,t)f(y,t)f(y,t) 即 y 在 t 处的导数以函数的形式表示y0: array 初始条件 y0y_0y0对于常微分方程组 y0y_0y0 则为数组向量t: array 求解函数值对应的时间点的序列。序列的第一个元素是与初始条件 y0y_0y0 对应的初始时间 t0t_0t0时间序列必须是单调递增或单调递减的允许重复值。
其它参数简介如下 args: 向导数函数 func 传递参数。当导数函数 f(y,t,p1,p2,..)f(y,t,p1,p2,..)f(y,t,p1,p2,..) 包括可变参数 p1,p2… 时通过 args (p1,p2,…) 可以将参数p1,p2… 传递给导数函数 func。argus 的用法参见 2.4 中的实例2。 Dfun: func 的雅可比矩阵行优先。如果 Dfun 未给出则算法自动推导。 col_deriv: 自动推导 Dfun的方式。 printmessg: 布尔值。控制是否打印收敛信息。 其它参数用于控制求解算法的参数一般情况可以忽略。
odeint 的主要返回值
y: array 数组形状为 (len(t),len(y0)给出时间序列 t 中每个时刻的 y 值。 5. 求洛伦兹Lorenz方程的数值解
6.1 例题 2求洛伦兹Lorenz方程的数值解
洛伦兹Lorenz混沌吸引子的轨迹可以由如下的 3个微分方程描述
{dxdtσ(y−x)dydtx(ρ−z)−ydzdtxy−βz\begin{cases} \begin{aligned} \frac{dx}{dt} \sigma (y-x)\\ \frac{dy}{dt} x (\rho-z) - y\\ \frac{dz}{dt} xy - \beta z\\ \end{aligned} \end{cases} ⎩⎨⎧dtdxσ(y−x)dtdyx(ρ−z)−ydtdzxy−βz
洛伦兹方程将大气流体运动的强度 x 与水平和垂直方向的温度变化 y 和 z 联系起来进行大气对流系统的模拟现已广泛应用于天气预报、空气污染和全球气候变化的研究。参数 σ\sigmaσ 称为普兰特数ρ\rhoρ 是规范化的瑞利数β\betaβ 和几何形状相关。洛伦兹方程是非线性微分方程组无法求出解析解只能使用数值方法求解。 6.2 洛伦兹Lorenz方程问题的编程步骤
以该题为例讲解 scipy.integrate.odeint() 求解常微分方程初值问题的步骤 导入 scipy、numpy、matplotlib 包 定义导数函数 lorenz(W, t, p, r, b) 注意 odeint() 函数中定义导数函数的标准形式是 f(y,t)f(y,t)f(y,t) 对于微分方程组 y 表示向量。 为避免混淆我们记为 W[x,y,z]W[x,y,z]W[x,y,z]函数 lorenz(W,t) 定义导数函数 f(W,t)f(W,t)f(W,t) 。 用 p,r,b 分别表示方程中的参数 σ、ρ、β\sigma、\rho、\betaσ、ρ、β则对导数定义函数编程如下
# 导数函数求 W[x,y,z] 点的导数 dW/dt
def lorenz(W,t,p,r,b):x, y, z W # W[x,y,z]dx_dt p*(y-x) # dx/dt p*(y-x), p: sigmady_dt x*(r-z) - y # dy/dt x*(r-z)-y, r:rhodz_dt x*y - b*z # dz/dt x*y - b*z, b;betareturn np.array([dx_dt,dy_dt,dz_dt])定义初值 W0W_0W0 和 WWW 的定义区间 [t0,t][t_0,\ t][t0, t] 调用 odeint() 求 WWW 在定义区间 [t0,t][t_0,\ t][t0, t] 的数值解。 注意例程中通过 argsparas 或 args (10.0,28.0,3.0) 将参数 (p,r,b) 传递给导数函数 lorenz(W,t,p,r,b)。参数 (p,r,b) 当然也可以不作为函数参数传递而是在导数函数 lorenz() 中直接设置。但例程的参数传递方法使导数函数结构清晰、更为通用。另外对于可变参数问题使用这种参数传递方式就非常方便。 6.3 洛伦兹Lorenz方程问题 Python 例程
# 2. 求解微分方程组初值问题(scipy.integrate.odeint)
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D# 导数函数, 求 W[x,y,z] 点的导数 dW/dt
def lorenz(W,t,p,r,b): # by youcansx, y, z W # W[x,y,z]dx_dt p*(y-x) # dx/dt p*(y-x), p: sigmady_dt x*(r-z) - y # dy/dt x*(r-z)-y, r:rhodz_dt x*y - b*z # dz/dt x*y - b*z, b;betareturn np.array([dx_dt,dy_dt,dz_dt])t np.arange(0, 30, 0.01) # 创建时间点 (start,stop,step)
paras (10.0, 28.0, 3.0) # 设置 Lorenz 方程中的参数 (p,r,b)# 调用ode对lorenz进行求解, 用两个不同的初始值 W1、W2 分别求解
W1 (0.0, 1.00, 0.0) # 定义初值为 W1
track1 odeint(lorenz, W1, t, args(10.0, 28.0, 3.0)) # args 设置导数函数的参数
W2 (0.0, 1.01, 0.0) # 定义初值为 W2
track2 odeint(lorenz, W2, t, argsparas) # 通过 paras 传递导数函数的参数# 绘图
fig plt.figure()
ax Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2], colormagenta) # 绘制轨迹 1
ax.plot(track2[:,0], track2[:,1], track2[:,2], colordeepskyblue) # 绘制轨迹 2
ax.set_title(Lorenz attractor by scipy.integrate.odeint)
plt.show()6.4 洛伦兹Lorenz方程问题 Python 例程运行结果 【本节完】
