数学建模:微分方程模型—常微分方程数值解算法及 Python 实现
目录
- 一、显式欧拉 (Euler) 法
- 二、显式欧拉法的改进
- 隐式欧拉法 (后退欧拉法)
- 梯形法
- 两步欧拉法 (中点法)
- 预报 - 校正法 (改进欧拉法)
- 三、龙格 - 库塔 (Runge-Kutta) 法
- 基本思路
- 2 阶龙格 - 库塔法
- 4 阶经典龙格 - 库塔法 (最常用)
- 一般的龙格 - 库塔法
考虑一阶常微分方程的初值问题
{dydx=f(x,y)x∈[a,b]y(a)=y0{\begin{cases} {\frac{{dy}}{{dx}} = f(x,y)\quad x \in [a,b]} \\ {y(a) = {y_0}} \end{cases}}{dxdy=f(x,y)x∈[a,b]y(a)=y0
只要 f(x,y)f(x, y)f(x,y) 在 [a,b]×R1[a, b] \times \mathbb{R}^{1}[a,b]×R1 上连续,且关于 yyy 满足 Lipschitz 条件,即存在与 x,yx, yx,y 无关的常数 LLL, s.t.s.t.s.t. ∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣\left|f\left(x, y_{1}\right)-f\left(x, y_{2}\right)\right| \leq L\left|y_{1}-y_{2}\right|∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣, 对任意定义在 [a,b][a, b][a,b] 上的 y1(x)y_{1}(x)y1(x) 和 y2(x)y_{2}(x)y2(x) 成立,则上述初值问题的连续可微解 y(x)y(x)y(x) 在 [a,b][a, b][a,b] 上存在且唯一。
常微分方程数值解要计算出解函数 y(x)y(x)y(x) 在一系列节点 a=x0<x1<…<xn=ba=x_{0}<x_{1}<\ldots<x_{n}=ba=x0<x1<…<xn=b 处的近似值 yi≈y(xi)(i=1,…,n)y_{i} \approx y\left(x_{i}\right) \ (i=1, \ldots, n)yi≈y(xi) (i=1,…,n),节点间距 hi=xi+1−xi(i=0,…,n−1)h_{i}=x_{i+1}-x_{i}(i=0, \ldots, n-1)hi=xi+1−xi(i=0,…,n−1) 称为步长,通常采用等距节点,即取 hi=hh_i = hhi=h (常数)。
一、显式欧拉 (Euler) 法
在 xyx yxy 平面上, 微分方程的解 y=y(x)y=y(x)y=y(x) 称为积分曲线. 积分曲线上一点 (x,y)(x, y)(x,y) 的切线斜率等于 函数 f(x,y)f(x, y)f(x,y) 的值.
从初始点 P0(x0,y0)P_{0}\left(x_{0}, y_{0}\right)P0(x0,y0) 出发, 沿切线方向推进到 x=x1x=x_{1}x=x1 上一点 P1P_{1}P1, 再从 P1P_{1}P1 沿切线方向推进到 x=x2x=x_{2}x=x2 上一点 P2P_{2}P2, 循此前进作出一条折线 P0P1P2⋯P_{0} P_{1} P_{2} \cdotsP0P1P2⋯, 作为积分曲线的一个近似. 即向前取差商近似导数
{y0=y(x0),yi+1=yi+hf(xi,yi)(i=0,1,⋯,n−1)\begin{cases} {y_0} = y({x_0}),\\ {y_{i + 1}} = {y_i} + {h}f({x_i},{y_i}) \quad (i = 0,1, \cdots ,n - 1) \end{cases} {y0=y(x0),yi+1=yi+hf(xi,yi)(i=0,1,⋯,n−1)
称为 (显式) 欧拉法
定义 在假设 yi=y(xi)y_{i}=y\left(x_{i}\right)yi=y(xi) ,即第 iii 步计算是精确的前提下,Ri=y(xi+1)−yi+1R_{i}=y\left(x_{i+1}\right)-y_{i+1}Ri=y(xi+1)−yi+1 称为局部截断误差。
定义 若某算法的局赔截断误差为 O(hp+1)O\left(h^{p+1}\right)O(hp+1) ,则称该算法有 ppp 阶精度。
显式欧拉格式具有 1 阶精度
Python 实现:
import numpy as np
#显式欧拉法的Python实现
def Euler(f,x0,y0,h,l):#f函数, x0,y0初值, h步长, l数量xn, yn = x0, y0n = 0ns, xs, ys = [n], [x0], [y0]while n <= l:n += 1xn += hyn = yn + f(xn,yn) * hns.append(n)xs.append(xn)ys.append(yn)return (ns,xs,ys)
二、显式欧拉法的改进
隐式欧拉法 (后退欧拉法)
将常微分方程化为积分方程
y(xi+1)=y(xi)+∫xixi+1f(x,y(x))dx.{y}\left(x_{i+1}\right)= {y}\left(x_{i}\right)+\int_{x_{i}}^{x_{i+1}} {f}(x, {y}(x)) {d} x . y(xi+1)=y(xi)+∫xixi+1f(x,y(x))dx.
采用左矩形公式计算积分
yi+1=yi+hf(xi,yi),{y}_{i+1}= {y}_{i}+h {f}\left(x_{i}, {y}_{i}\right), yi+1=yi+hf(xi,yi),
即是显式欧拉法的递推公式.
用右矩形公式计算积分, 即是向后取差商近似导数
yi+1=yi+hf(xi+1,yi+1)(i=0,…,n−1)y_{i+1}=y_{i}+h f\left(x_{i+1}, y_{i+1}\right) \quad(i=0, \ldots, n-1) yi+1=yi+hf(xi+1,yi+1)(i=0,…,n−1)
称为隐式欧拉法或后退欧拉法.
隐式欧拉格式具有 1 阶精度
梯形法
即显、隐式两种算法的平均
yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+1)](i=0,...,n−1){y_{i + 1}} = {y_i} + \frac{h}{2}[f({x_i},{y_i}) + f({x_{i + 1}},{y_{i + 1}})]\quad (i = 0,\;...\;,n - 1) yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+1)](i=0,...,n−1)
称为梯形法.
梯形格式具有 2 阶精度
两步欧拉法 (中点法)
取中心差商近似导数
yi+1=yi−1+2hf(xi,yi)(i=1,...,n−1){y_{i + 1}} = {y_{i - 1}} + 2h\,f({x_i},{y_i})\quad (i = 1,\;...\;,n - 1) yi+1=yi−1+2hf(xi,yi)(i=1,...,n−1)
称为两步欧拉法或中点法.
中点格式具有 2 阶精度
预报 - 校正法 (改进欧拉法)
先用显式欧拉法作预报, 算出 yˉi+1=yi+hf(xi,yi)\bar{y}_{i+1}=y_{i}+h f\left(x_{i}, y_{i}\right)yˉi+1=yi+hf(xi,yi)
再将 yˉi+1\bar{y}_{i+1}yˉi+1 代入隐式梯形法的右边作校正,得到
yi+1=yi+h2[f(xi,yi)+f(xi+1,yˉi+1)]y_{i+1}=y_{i}+\frac{h}{2}\left[f\left(x_{i}, y_{i}\right)+f\left(x_{i+1}, \bar{y}_{i+1}\right)\right] yi+1=yi+2h[f(xi,yi)+f(xi+1,yˉi+1)]
即
yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+hf(xi,yi))](i=0,…,n−1)y_{i+1}=y_{i}+\frac{h}{2}\left[f\left(x_{i}, y_{i}\right)+f\left(x_{i+1}, y_{i}+h f\left(x_{i}, y_{i}\right)\right)\right] \quad(i=0, \ldots, n-1) yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+hf(xi,yi))](i=0,…,n−1)
预报 - 校正法通常只迭代一两次就转人下一步的计算, 避免了梯形法和中点法每次迭代, 都要重新计算函数 f(x,y){f}(x, {y})f(x,y), 减少了计算量.
预报 - 校正法具有 2 阶精度
三、龙格 - 库塔 (Runge-Kutta) 法
基本思路
积分方程右端积分的求积公式可以表示为
∫xixi+1f(x,y(x))dx≈h∑i=1rcif(xi+λih,y(xi+λih))\int_{x_{i}}^{x_{i+1}} {f}(x, {y}(x)) \mathrm{d} x \approx h \sum_{i=1}^{r} c_{i} {f}\left(x_{i}+\lambda_{i} h, {y}\left(x_{i}+\lambda_{i} h\right)\right) ∫xixi+1f(x,y(x))dx≈hi=1∑rcif(xi+λih,y(xi+λih))
将上式表示为
yn+1=yn+hϕ(xn,yn,h){y}_{n+1}= {y}_{n}+h {\phi}\left(x_{n}, {y}_{n}, h\right) yn+1=yn+hϕ(xn,yn,h)
其中
ϕ(xi,yi,h)=∑i=1rciKiK1=f(xi,yi)Ki=f(xi+λih,yi+h∑j=1i−1μijKj)(i=2,⋯,r)\begin{gathered} &\phi\left(x_{i}, {y}_{i}, h\right)=\sum_{i=1}^{r} c_{i} {K}_{i} \\ & {K}_{1}= {f}\left(x_{i}, {y}_{i}\right) \\ & {K}_{i}= {f}\left(x_{i}+\lambda_{i} h, {y}_{i}+h \sum_{j=1}^{i-1} \mu_{i j} {K}_{j}\right) \quad(i=2, \cdots, r) \end{gathered} ϕ(xi,yi,h)=i=1∑rciKiK1=f(xi,yi)Ki=f(xi+λih,yi+hj=1∑i−1μijKj)(i=2,⋯,r)
称为 rrr 阶显式龙格 - 库塔法. 当 r=1r=1r=1 时, ϕ(xi,yi,h)=f(xi,yi)\phi\left(x_{i}, y_{i}, h\right)=f\left(x_{i}, y_{i}\right)ϕ(xi,yi,h)=f(xi,yi), 就是欧拉法.
2 阶龙格 - 库塔法
取 r=2r=2r=2 时
{yi+1=yi+h[λ1K1+λ2K2]K1=f(xi,yi)K2=f(xi+ph,yi+phK1)\left\{\begin{array}{l} y_{i+1}=y_{i}+h\left[\lambda_{1} K_{1}+\lambda_{2} K_{2}\right] \\ K_{1}=f\left(x_{i}, y_{i}\right) \\ K_{2}=f\left(x_{i}+p h, y_{i}+p h K_{1}\right) \end{array}\right. ⎩⎨⎧yi+1=yi+h[λ1K1+λ2K2]K1=f(xi,yi)K2=f(xi+ph,yi+phK1)
其中 λ1+λ2=1,λ2p=12{\lambda _1} + {\lambda_2} = 1, \ {\lambda_2}p = \frac{1}{2}λ1+λ2=1, λ2p=21, 算法有 2 阶精度
若取 λ1=0,λ2=1{\lambda _1}=0, {\lambda_2} = 1λ1=0,λ2=1, 则 p=12p = \frac{1}{2}p=21, 则
{yi+1=yi+hK2K1=f(xi,yi)K2=f(xi+h2,yi+h2K1)\left\{\begin{array}{l} {y}_{i+1}= {y}_{i}+h {K}_{2} \\ {K}_{1}= {f}\left(x_{i}, {y}_{i}\right) \\ {K}_{2}= {f}\left(x_{i}+\frac{h}{2}, {y}_{i}+\frac{h}{2} {K}_{1}\right) \end{array}\right. ⎩⎨⎧yi+1=yi+hK2K1=f(xi,yi)K2=f(xi+2h,yi+2hK1)
即是中点法
4 阶经典龙格 - 库塔法 (最常用)
4 阶龙格 - 库塔法是最常用的龙格 - 库塔法
{yi+1=yi+h6(K1+2K2+2K3+K4)K1=f(xi,yi)K2=f(xi+h2,yi+h2K1)K3=f(xi+h2,yi+h2K2)K4=f(xi+h,yi+hK3)\left\{\begin{array}{l} y_{i+1}=y_{i}+\frac{h}{6}\left(K_{1}+2 K_{2}+2 K_{3}+K_{4}\right) \\ K_{1}=f\left(x_{i}, y_{i}\right) \\ K_{2}=f\left(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2} K_{1}\right) \\ K_{3}=f\left(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2} K_{2}\right) \\ K_{4}=f\left(x_{i}+h, y_{i}+h K_{3}\right) \end{array}\right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧yi+1=yi+6h(K1+2K2+2K3+K4)K1=f(xi,yi)K2=f(xi+2h,yi+2hK1)K3=f(xi+2h,yi+2hK2)K4=f(xi+h,yi+hK3)
Python 实现:
import numpy as np
#四阶龙格库塔法的Python实现
def RK4(f,x0,y0,h,maxx):#f函数, x0,y0初值, h步长, maxx: x最大值xn, yn = x0, y0n = 0ns, xs, ys = [n], [xn], [yn]while xn < maxx:n += 1xn += hk1 = f(xn,yn)k2 = f(xn + (h / 2),yn + (h * k1) / 2)k3 = f(xn + (h / 2),yn + (h * k2) / 2)k4 = f(xn + h,yn + h * k3)yn = yn + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)ns.append(n)xs.append(xn)ys.append(yn)return (ns,xs,ys)
一般的龙格 - 库塔法
{yi+1=yi+h[λ1K1+λ2K2+…+λmKm]K1=f(xi,yi)K2=f(xi+α2h,yi+β21hK1)K3=f(xi+α3h,yi+β31hK1+β32hK2)⋯⋯Km=f(xi+αmh,y+βm1hK1+βm2hK2+…+βmm−1hKm−1)\left\{\begin{aligned} y_{i+1} &=y_{i}+h\left[\lambda_{1} K_{1}+\lambda_{2} K_{2}+\ldots+\lambda_{m} K_{m}\right] \\ K_{1} &=f\left(x_{i}, y_{i}\right) \\ K_{2} &=f\left(x_{i}+\alpha_{2} h, y_{i}+\beta_{21} h K_{1}\right) \\ K_{3} &=f\left(x_{i}+\alpha_{3} h, y_{i}+\beta_{31} h K_{1}+\beta_{32} h K_{2}\right) \\ & \cdots \cdots \\ K_{m} &=f\left(x_{i}+\alpha_{m} h, y+\beta_{m 1} h K_{1}+\beta_{m 2} h K_{2}+\ldots+\beta_{m m-1} h K_{m-1}\right) \end{aligned}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧yi+1K1K2K3Km=yi+h[λ1K1+λ2K2+…+λmKm]=f(xi,yi)=f(xi+α2h,yi+β21hK1)=f(xi+α3h,yi+β31hK1+β32hK2)⋯⋯=f(xi+αmh,y+βm1hK1+βm2hK2+…+βmm−1hKm−1)
由于龙格-库塔法的导出基于泰勒展开, 故精度主要受解函数的光滑性影响. 对于光滑性不太好的解, 最好采用低阶算法而将步长 hhh 取小.
数学建模:微分方程模型—常微分方程数值解算法及 Python 实现相关推荐
- 数学建模——微分方程模型的求解
文章目录 微分方程的符号解法 微分方程数值解法 一些常用的微分方程模型(学习中,持续更新) Logistics模型 传染病模型 本文介绍微分方程的求解,不介绍微分方程的建立方法 微分方程的符号解法 求 ...
- Python小白的数学建模课-17.条件最短路径算法
条件最短路径问题,指带有约束条件.限制条件的最短路径问题.例如: 顶点约束,包括必经点或禁止点的限制: 边的约束,包括必经路段.禁行路段和单向路段:无权路径长度的限制,如要求经过几步或不超过几步到达终 ...
- 数学建模常用模型(一):灰色预测法
数学建模常用模型(一):灰色预测法 灰色预测法是一种用于处理少量数据.数据质量较差或者缺乏历史数据的预测方法.它适用于一些非线性.非平稳的系统,尤其在短期预测和趋势分析方面有着广泛的应用.灰色预测法作 ...
- 数学建模常见模型总结
数学建模常见模型总结 一.插值 当已有数据量不够,需要补充,且认定已有数据可信时,通常利用函数插值方法. 常用插值方法 拉格朗日插值 分段线性插值 Hermite 三次样条插值 克里金法 matlab ...
- 数学建模——支持向量机模型详解Python代码
数学建模--支持向量机模型详解Python代码 from numpy import * import random import matplotlib.pyplot as plt import num ...
- 数学建模——线性规划模型详解Python代码
数学建模--线性规划模型详解Python代码 标准形式为: min z=2X1+3X2+x s.t x1+4x2+2x3>=8 3x1+2x2>=6 x1,x2,x3>=0 上述线性 ...
- 数学建模传染病模型_数学建模| 时间序列模型
1 数学建模 时间序列模型 1.与实践有关系的一组数据,叫做时间序列: 2.得到时间序列的数据后,要构建模型,其中平稳时间序列的模型,是本节课重点介绍的: 3.y=at+季节性+周期性 一. ...
- 数学建模常用模型04:灰色关联分析法
数学建模常用模型04:灰色关联分析法 灰色关联分析法 本文所用的资料参考来源:美赛资料网:美赛资料网 与灰色预测模型一样,比赛不能优先使用,灰色关联往往可以与层次分析结合使用.层次分析用在确定权重上面 ...
- 华中杯 数学建模 A题简单复盘(附Python源码)
华中杯 A题简单复盘(附Python 源码) 文章目录 华中杯 A题简单复盘(附Python 源码) 前言 题目简介 问题背景 题目以及思路 分批算法设计 MindMap 遗传算法优缺点 优点 缺点 ...
最新文章
- 挥手送别 2019,翘首期待 2020
- 用gcc gvim编译程序
- Python rjust() 方法
- 低学历程序员的红利来了,这个政策来的太惊喜!
- mac 硬盘未推出 硬盘无法读取_在Mac上(正确的)格式化U盘
- 公布自己的pods到CocoaPods trunk 及问题记录
- 函数参数的值传递和地址传递
- 加载文件流_jvm类加载的过程
- Android之创建选项菜单
- GaussDB(openGauss)宣布开源,性能超越 MySQL 与 PostgreSQL
- ASPxGridView1用法-
- Facebook又开两处AI实验室,在西雅图和匹兹堡招兵买马
- element tree不刷新视图_安卓从入门到进阶第五章(视图查看)
- linux下常用vim命令
- win10安装VS2008失败解决方案
- 无线专题 路由器和交换机、光猫的区别
- 国外、国内Hadoop的应用现状
- 3dmax 2014加载panda3d插件失败
- JAVA12_10总结
- java高效随机生成随机(英文+数字),可自定义
热门文章
- 阿里云轻量应用服务器开启minecraft基岩版服务器(bedrock)
- vue项目对接pad端——混合开发总结
- 研究表明:菜鸟爱用右脑,专家爱用左脑!
- table html 合并列,html table上下行合并
- 基于北斗导航定位系统的设计与实现(论文+程序设计源码+数据库文件)
- 多通路fpga 通信_【论文精选】基于FPGA的EtherCAT从站通信链路分析与验证
- python中的高级特性
- spark TF-IDF特征提取生成文章关键词
- 如何在微信复制链接直接可以用浏览器打开 微信调用手机浏览器打开指定链接
- 祝你一路顺风_吴奇隆_酷音小伟编曲_C调简单版