当前位置:   article > 正文

数学建模:微分方程模型— Python 求解_python欧拉法求解微分方程

python欧拉法求解微分方程

使用 Python 求常微分方程的数值求解通常是基于一阶方程进行的,高阶微分方程要化成一阶方程组。

例:使用显式欧拉法和四阶龙格库塔法计算Lorenz模型

Lorenz模型是由美国气象学家Lorenz在研究大气运动时,通过简化对流模型,只保留3个变量提出的一个完全确定性的一阶自治常微分方程组(不显含时间变量),其方程为

{ x ˙ = σ ( y − x ) , y ˙ = ρ x − y − x z , z ˙ = x y − β z .

{x˙=σ(yx),y˙=ρxyxz,z˙=xyβz.
x˙=σ(yx),y˙=ρxyxz,z˙=xyβz.

其中,参数 σ \sigma σ 为 Prandtl 数, ρ \rho ρ 为 Rayleigh 数, β \beta β 为方向比。

这里分别用显式欧拉法和四阶龙格库塔法给出 Lorenz 模型在 σ = 10 , ρ = 28 , β = 8 3 \sigma=10, \rho=28, \beta=\frac{8}{3} σ=10,ρ=28,β=38 时系统的三维演化轨迹,和系统从两个靠的很近的初值出发(相差仅 0.0001 0.0001 0.0001)后,解的偏差演化曲线。显示欧拉法及四阶龙格库塔法详见此处

#显示欧拉法解lorenz方程
import numpy as np
import matplotlib.pyplot as plt

#显式欧拉法
def Euler(f,x0,y0,h,l):
    #f函数, x0,y0初值, h步长, l数量
    xn, yn = x0, y0
    n = 0
    ys = [y0]
    while n <= l:
        n += 1
        xn += h
        yn = yn + f(xn,yn) * h
        ys.append(yn)
    return np.array(ys)

#定义函数
def lorenz(t,w):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
#t=np.arange(0, 50, 0.01)  #创建时间点
sol1=Euler(lorenz, 0, [0.0, 1.0, 0.0], 0.01, 5000)  #第一个初值问题求解
sol2=Euler(lorenz, 0, [0.0, 1.0001, 0.0], 0.01, 5000)  #第二个初值问题求解

#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1[:,0], sol1[:,1], sol1[:,2],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1[:,0]-sol2[:,0], sol1[:,1]-sol2[:,1], sol1[:,2]-sol2[:,2],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1, '\n\n', "sol1-sol2=", sol1-sol2)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
'
运行

欧拉l

import numpy as np
import matplotlib.pyplot as plt

#四阶龙格库塔
def RK4(f,x0,y0,h,maxx):
    #f函数, x0,y0初值, h步长, maxx: x最大值
    xn, yn = x0, y0
    n = 0
    ys = [yn]
    while xn < maxx:
        n += 1
        xn += h
        k1 = 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)
        ys.append(yn)
    return np.array(ys)

#定义函数
def lorenz(t,w):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
sol1=RK4(lorenz, 0, [0.0, 1.0, 0.0], 0.01, 50)  #第一个初值问题求解
sol2=RK4(lorenz, 0, [0.0, 1.0001, 0.0], 0.01, 50)  #第二个初值问题求解

#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1[:,0], sol1[:,1], sol1[:,2],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1[:,0]-sol2[:,0], sol1[:,1]-sol2[:,1], sol1[:,2]-sol2[:,2],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1, '\n\n', "sol1-sol2=", sol1-sol2)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
'
运行

KLL
两个程序的输出图像中,左图显示出我们经常听到的“蝴蝶效应”。右图给出了系统从两个靠的很近的初值出发(相差仅0.0001)后,解的偏差演化曲线,从图中看出随着时间的增大,两个解的差异越来越大,这正是动力系统对初值敏感性的直观表现。

scipy.integrate.odeint 求解微分方程模型

scipy.integrate 模块的 odeint 函数求常微分方程的数值解的基本调用格式为:

from scipy.integrate import odeint
sol=odeint(func, y0, t)
  • 1
  • 2

其中func是定义微分方程的函数或匿名函数,y0是初始条件的序列,t是一个自变量取值的序列(t的第一个元素一定为初始时刻),返回值sol是对应于序列t中元素的数值解,如果微分方程组中有 n n n个函数,返回值sol是 n n n列的矩阵,第 i ( i = 1 , 2 , ⋯   , n ) i(i = 1,2, \cdots ,n) i(i=1,2,,n)列对应于第 i i i个函数的数值解。

示例程序: 用 scipy.integrate.odeint 求解 Lorenz 方程

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
#定义函数
def lorenz(w,t):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
t=np.arange(0, 50, 0.01)  #创建时间点
sol1=odeint(lorenz, [0.0, 1.0, 0.0], t)  #第一个初值问题求解
sol2=odeint(lorenz, [0.0, 1.0001, 0.0], t)  #第二个初值问题求解
#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1[:,0], sol1[:,1], sol1[:,2],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1[:,0]-sol2[:,0], sol1[:,1]-sol2[:,1], sol1[:,2]-sol2[:,2],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1, '\n\n', "sol1-sol2=", sol1-sol2)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23

scipy.integrate.solve_ivp 求解微分方程模型

求解 { d y d t = f ( t , y ) , y ( t 0 ) = y 0.

{dydt=f(t,y),y(t0)=y0.
dtdy=f(t,y),y(t0)=y0.

调用方式:

from scipy.integrate import solve_ivp
slove=scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, 
		dense_output=False, events=None, vectorized=False, args=None, **options)
  • 1
  • 2
  • 3

参数:

  1. fun: 定义微分方程的函数, fun(t,y)
  2. t_span: 2 元组的浮点数, 积分区间 (t0, tf)
  3. y0: 数组, 形状 (n,), 初始状态
  4. method: 使用的积分方法, 如Runge-Kutta方法 (‘RK23’、‘RK45’、‘DOP853’)
  5. t_eval: 数组 或 None, 可选, 求解的时间, 必须在 t_span 内, 如果无(默认),则使用求解器选择的点

返回值:

print(slove.t) #时间
print{slove.y) #解
  • 1
  • 2

示例程序: 用 scipy.integrate.odeint 求解 Lorenz 方程

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
#定义函数
def lorenz(t,w):
    sigma=10; rho=28; beta=8/3
    x, y, z=w
    return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
t=np.arange(0, 50, 0.01)  #创建时间点
sol1=solve_ivp(lorenz, (0,50), np.array([0.0, 1.0, 0.0]), 'RK45', t)  #第一个初值问题求解
sol2=solve_ivp(lorenz, (0,50), np.array([0.0, 1.0001, 0.0]), 'RK45', t)  #第二个初值问题求解
#绘图
plt.rc('font',size=16); plt.rc('text',usetex=True)
#图A
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1.y[0,:], sol1.y[1,:], sol1.y[2,:],'r')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
#图B
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1.y[0,:]-sol2.y[0,:], sol1.y[1,:]-sol2.y[1,:], sol1.y[2,:]-sol2.y[2,:],'g')
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$'); ax2.set_zlabel('$z$')
plt.show()
#print("sol1=",sol1.y, '\n\n', "sol1-sol2=", sol1.y-sol2.y)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/正经夜光杯/article/detail/814108
推荐阅读
相关标签
  

闽ICP备14008679号