python如何求解微分方程_常微分方程数值解:Python求解
这里对使用python求解常微分方程提供两种思路,一种是自己编程实现欧拉法,改进欧拉法或者四阶龙格库塔,这样有助于理解上述三种数值计算方法的原理;一种是调用python已有的库,不再重复造轮子。
本文对上述两种思路都给出代码示例,并进行比较;同时针对单个微分方程和含有多个微分方程的微分方程组给出代码示例
1. 常微分方程定义
凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程。
未知函数是一元函数的微分方程称作常微分方程
未知数是多元函数的微分方程称作偏微分方程。
微分方程中出现的未知函数最高阶导数的阶数,称为微分方程的阶数。
2. 调用现有的库
scipy中提供了用于解常微分方程的函数odeint(),完整的调用形式如下:
scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, \
rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0,hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
实际使用中,还是主要使用前三个参数,即微分方程的描写函数、初值和需要求解函数值对应的的时间点。接收数组形式。这个函数,要求微分方程必须化为标准形式,即
,实际操作中会发现,高阶方程的标准化工作,其实是解微分方程最主要的工作。
示例1:单个微分方程
import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def func(y, t):
return t * math.sqrt(y)
YS=odeint(func,y0=1,t=np.arange(0,10.1,0.1))
t=np.arange(0,10.1,0.1)
plt.plot(t, YS, label='odeint')
plt.legend()
plt.show()
结果如下:
图.1
示例2:微分方程组
与单个微分方程不同的是,此时的函数变成了向量函数
from scipy.integrate import odeint
import numpy as np
def lorenz(w, t, p, r, b):
# 给出位置矢量w,和三个参数p, r, b计算出
# dx/dt, dy/dt, dz/dt的值
x, y, z = w
# 直接与lorenz的计算公式对应
return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
t = np.arange(0, 30, 0.01) # 创建时间点
# 调用ode对lorenz进行求解, 用两个不同的初始值
track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))
# 绘图
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2])
ax.plot(track2[:,0], track2[:,1], track2[:,2])
plt.show()
结果如下
图.2
3. 自己编程实现欧拉法/改进欧拉法/四阶龙格库塔
示例 3:单个函数使用四阶龙格库塔
import math
import numpy as np
import matplotlib.pyplot as plt
def runge_kutta(y, x, dx, f):
""" y is the initial value for y
x is the initial value for x
dx is the time step in x
f is derivative of function y(t)
"""
k1 = dx * f(y, x)
k2 = dx * f(y + 0.5 * k1, x + 0.5 * dx)
k3 = dx * f(y + 0.5 * k2, x + 0.5 * dx)
k4 = dx * f(y + k3, x + dx)
return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.
if __name__=='__main__':
t = 0.
y = 1.
dt = .1
ys, ts = [], []
def func(y, t):
return t * math.sqrt(y)
while t <= 10:
y = runge_kutta(y, t, dt, func)
t += dt
ys.append(y)
ts.append(t)
exact = [(t ** 2 + 4) ** 2 / 16. for t in ts]
plt.plot(ts, ys, label='runge_kutta')
plt.plot(ts, exact, label='exact')
plt.legend()
#plt.show()
结果如下:
图.3
示例4:示例1和示例3放在一起进行对比
import math
import numpy as np
import matplotlib.pyplot as plt
def runge_kutta(y, x, dx, f):
""" y is the initial value for y
x is the initial value for x
dx is the time step in x
f is derivative of function y(t)
"""
k1 = dx * f(y, x)
k2 = dx * f(y + 0.5 * k1, x + 0.5 * dx)
k3 = dx * f(y + 0.5 * k2, x + 0.5 * dx)
k4 = dx * f(y + k3, x + dx)
return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.
if __name__=='__main__':
t = 0.
y = 1.
dt = .1
ys, ts = [], []
def func(y, t):
return t * math.sqrt(y)
while t <= 10:
y = runge_kutta(y, t, dt, func)
t += dt
ys.append(y)
ts.append(t)
from scipy.integrate import odeint
YS=odeint(func,y0=1, t=np.arange(0,10.1,0.1))
plt.plot(ts, ys, label='runge_kutta')
plt.plot(ts, YS, label='odeint')
plt.legend()
plt.show()
结果如下
图.4
示例5:多个微分方程(欧拉法)
import numpy as np
"""
移动方程:
t时刻的位置P(x,y,z)
steps:dt的大小
sets:相关参数
"""
def move(P, steps, sets):
x, y, z = P
sgima, rho, beta = sets
# 各方向的速度近似
dx = sgima * (y - x)
dy = x * (rho - z) - y
dz = x * y - beta * z
return [x + dx * steps, y + dy * steps, z + dz * steps]
# 设置sets参数
sets = [10., 28., 3.]
t = np.arange(0, 30, 0.01)
# 位置1:
P0 = [0., 1., 0.]
P = P0
d = []
for v in t:
P = move(P, 0.01, sets)
d.append(P)
dnp = np.array(d)
# 位置2:
P02 = [0., 1.01, 0.]
P = P02
d = []
for v in t:
P = move(P, 0.01, sets)
d.append(P)
dnp2 = np.array(d)
"""
画图
"""
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(dnp[:, 0], dnp[:, 1], dnp[:, 2])
ax.plot(dnp2[:, 0], dnp2[:, 1], dnp2[:, 2])
plt.show()
结果如下:
图.5
参考:
python如何求解微分方程_常微分方程数值解:Python求解相关推荐
- 数学建模:微分方程模型—常微分方程数值解算法及 Python 实现
目录 一.显式欧拉 (Euler) 法 二.显式欧拉法的改进 隐式欧拉法 (后退欧拉法) 梯形法 两步欧拉法 (中点法) 预报 - 校正法 (改进欧拉法) 三.龙格 - 库塔 (Runge-Kutta ...
- python如何求解微分方程_用Python数值求解偏微分方程
1 引言 微分方程是描述一个系统的状态随时间和空间演化的最基本的数学工具之一,其在物理.经济.工程.社会等各方面都有及其重要的应用.然而,只有很少的微分方程可以解析求解,尤其对于偏微分方程,能解析求解 ...
- 四阶龙格库塔法求解一次常微分方程组(python实现)
四阶龙格库塔法求解一次常微分方程组 一.前言 二.RK4求解方程组的要点 1. 将方程组转化为RK4求解要求的标准形式 2. 注意区分每个方程的独立性 三.python实现RK4求解一次常微分方程组 ...
- c++调用cplex求解例子_视频教程 | 用Python玩转运筹优化求解器IBM CPLEX(二)
编者按 优化求解器对于做运筹学应用的学生来说,意义重大. 然而直到今天,放眼望去,全网(包括墙外)几乎没有一个系统的Cplex中文求解器教程. 作为华人运筹学的最大的社区,『运筹OR帷幄』 责无旁贷, ...
- 常微分方程数值解matlab欧拉,MATLAB实验报告_常微分方程数值解
manlab软件应用试验题目 专业 序号 姓名 日期 实验3 常微分方程数值解 [实验目的] 1.掌握用MATLAB求微分方程初值问题数值解的方法: 2.通过实例学习微分方程模型解决简化的实际问题: ...
- 缩进对于python程序至关重要吗_缩进对于Python程序至关重要。
[判断题]用截面法求内力时,可以保留截开后构件的任意部分进行平衡计算 . [单选题]对大量生产的零件,机加工应采用( ) [单选题]对事件重要与紧急情况的分析,也是对: [单选题]在下列传动中,平均传 ...
- python变量定义大全_详解python变量与数据类型
这篇文章我们学习 Python 变量与数据类型 变量 变量来源于数学,是计算机语言中能储存计算结果或能表示值抽象概念,变量可以通过变量名访问.在 Python 中 变量命名规定,必须是大小写英文,数字 ...
- python金融量化书籍_超强干货 | Python金融数据量化分析教程+机器学习电子书
如今Python语言的学习已经上升到了国家战略的层面上.Python语言是人工智能的基础语言,国家相关教育部门对于"人工智能普及"格外重视,不仅将Python列入到小学.中学和高中 ...
- python积木式编程_实例讲解python函数式编程
函数式编程是使用一系列函数去解决问题,按照一般编程思维,面对问题时我们的思考方式是"怎么干",而函数函数式编程的思考方式是我要"干什么". 至于函数式编程的特点 ...
最新文章
- javascript语法糖_语法糖和JavaScript糖尿病
- Foundation ActionScript 3.0.With Flash CS3 And Flex ..
- 只想安安静静的做个程序员
- 用vim的方式操作你的软件
- 暖通空调系统计量表选型与应用
- 为什么阿里巴巴RPC接口不允许使用枚举类型?
- mockito 静态方法_Mockito –带有注释和静态方法的额外接口
- 【渝粤教育】电大中专金融与税收作业 题库
- 图片阴影怎么设置_HTML5 给图形绘制阴影
- 快速排序法实战入门(推荐)
- PHP操作tcpdf插件生成PDF
- Mac WebStorm 破解
- VMware 10激活码
- 最新第一波:全国信息化工程师软考-系统集成项目管理工程师(高级案例高分论文)
- 配置git的合并工具mergetool不生成.orig文件
- 获取公众号关注页面链接
- (回溯_04)组合总和
- Android各个版本代号及其特性
- 5G标准正式出炉 5G的杀手锏业务又在哪里呢
- 科技论文成功发表的技巧 (Get Published! Successful Scientific Writing)读书笔记