[Python知识库] python解常微分方程（组）

开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> Python知识库 -> python解常微分方程（组） -> 正文阅读

[Python知识库]python解常微分方程（组）

本人目前初三，能力所限，如有不足之处，还望多多指教。

一周前看到了一个视频，于是我便想用python来求解这个问题。

〇、分析

假设在平面内有一带电粒子q，质量为m。空间内存在匀强磁场B，方向垂直于平面向内即沿z轴负半轴，以及一个沿y轴负半轴的重力场。带电粒子从磁场内O点释放。

则可直接列出粒子的运动方程

$\large m\ddot{\vec{r}} = q \cdot \vec{v} \times \vec{B}$

将这个方程分解成x和y两个方向

$\large \begin{cases} \ddot{x} = -\frac{qB}{m}\dot{y}\\ \ddot{y} = -g+\frac{qB}{m}\dot{x}\\ \end{cases}$

联立即可求得该方程组的解。

一、sympy中的dsolve方法

#导入
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
init_printing()

首先声明符号x,y,q,m,B,g

#参量
q,m,B,t,g = symbols('q m B t g',real=True,nonzero=True,positive=True)
#函数
x,y = symbols('x y',real=True,nonzero=True,positive=True,cls=Function)

再将微分方程表示出来

eq1 = Eq(x(t).diff(t,t),-q*B*y(t).diff(t)/m)
eq2 = Eq(y(t).diff(t,t),-g+q*B*x(t).diff(t)/m)
sol = dsolve([eq1,eq2])

现在打印出sol(用jupyter notebook效果更好)?

sol

很明显，这个式子非常冗杂，用trigsimp()方法化简

x = trigsimp(sol[0].rhs)
y = trigsimp(sol[1].rhs)
x,y

然后，可以将里面的积分常数算出

#定义积分变量，避免报错
var('C1 C2 C3 C4')
#x(0)=0
e1 = Eq(x.subs({t:0}),0)
#x'(0)=0,要将subs放在diff后面
e2 = Eq(x.diff(t).subs({t:0}),0)
#y(0)=0
e3 = Eq(y.subs({t:0}),0)
#y'(0)=0
e4 = Eq(y.diff(t).subs({t:0}),0)
l = solve([e1,e2,e3,e4],[C1,C2,C3,C4])
l

紧接着，将积分常量替换到x和y里面,我们就得到了解的最终形式

x = x.subs(l)
y = y.subs(l)
x,y

当然，这个解还可以写成这种形式

$\large \begin{cases} x = \frac{gm^2}{B^2q^2}\left ( \frac{Bq}{m}t-sin\left ( \frac{Bq}{m}t \right ) \right )\\ y = \frac{gm^2}{B^2q^2}\left ( cos\left ( \frac{Bq}{m}t\right ) -1\right ) \end{cases}$

用plt画图来看看

ts = np.linspace(0,15,1000)
consts = {q:1,B:1,g:9.8,m:1}
fx = lambdify(t,x.subs(consts),'numpy')
fy = lambdify(t,y.subs(consts),'numpy')
plt.plot(fx(ts),fy(ts))
plt.grid()
plt.show()

但是sympy有一个缺点，当微分方程很复杂时，它会直接罢工。

于是另一个新的方法替我们解决了这个问题

二、scipy.integrate中的odeint方法

#导入
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from math import e

我们首先要创建一个函数，使它可以表示这个微分方程。

对于这个函数应该具备怎样的形式，先从一阶微分方程开始

比如说，我们要求解下面这个方程

$\large \frac{\mathrm{d}y }{\mathrm{d} x}+\frac{y}{x} = 0$

它的解通过一些简单的分离变量即得

$\large y = \frac{C}{x}$

下面用odeint来解出它，令y(1)=1

#一阶
def f(x,y):
  #将一阶导数用y和x构成的函数来表示
  dydx = -y/x
  return dydx
#初始条件
y0 = 1
#初值条件是在自变量x的下限取到。如y(1)就要生成一个从1开始的范围
x = np.linspace(1,5,100)
plt.xlim(0,5)
plt.ylim(0,5)
#odeint(）方法,tfirst属性指的是函数f的第一个参数是自变量
sol = odeint(f,y0,x,tfirst=True)
plt.plot(x,sol[:,0])
plt.grid()
plt.show()

同理，对于下面的一阶微分方程组

$\large \begin{cases} \frac{\mathrm{d}x }{\mathrm{d} t}=3x-2y\\ \frac{\mathrm{d}y }{\mathrm{d} t} = 2x - y + e^t \end{cases}$

此时我们设向量u=(x,y) （列向量），方程组化为

$\large \frac{\mathrm{d} \mathbf{u} }{\mathrm{d} x} = \binom{3x-2y}{2x-y+e^t}$

其实对于任意的一阶微分方程，都可以写成这样的形式

$\large \mathbf{u} = \left( \begin{matrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{matrix} \right ), f(t,\mathbf{u})= \left( \begin{matrix} f_1(t,x_1,x_2,\cdots,x_n) \\ f_2(t,x_1,x_2,\cdots,x_n) \\ \vdots \\ f_n(t,x_1,x_2,\cdots,x_n) \end{matrix} \right ), \frac{\mathrm{d} \mathbf{u}}{\mathrm{d} t} = \left( \begin{matrix} \frac{\mathrm{d}x_1 }{\mathrm{d} t} \\ \frac{\mathrm{d}x_2 }{\mathrm{d} t} \\ \vdots \\ \frac{\mathrm{d}x_n }{\mathrm{d} t} \end{matrix} \right )$

f的含义就是

$\large \\t \in (a,b)\\ \mathbf{u} \in \mathbf{D} \\ f:(a,b) \times \mathbf{D}\subseteq \mathbf{R} \times \mathbf{R^n} \rightarrow \mathbf{R^n}$

此时微分方程就可以化为

$\large \frac{\mathrm{d} \mathbf{u}}{\mathrm{d} t} = f(t,\mathbf{u})$

而f(t,u)正是我们要找的那个函数

odeint中支持向量输入,因此，可以这样构造这个函数

def f(t,u):
  #x1,x2...xn = u
  x,y = u
  #dx1dt=...,dx2dt=...,dxndt=...
  dxdt = 3*x-x*y
  dydt = 2*x-y+e**t
  #dudt = [dx1dt,dx2dt,...dxndt]
  dudt = [dxdt,dydt]
  return dudt
#x1(0)=...,x2(0)=...,xn(0)=...
u0=[0,0]
t = np.linspace(0,10,100)
sol = odeint(f,u0,t,tfirst=True)

因此，对于二阶常微分方程

$\large \frac{\mathrm{d^2}y }{\mathrm{d} x^2} + 4\frac{\mathrm{d}y }{\mathrm{d} x}+4y=e^x$

我们首先把这个方程化解成这样的形式

$\large \begin{cases} \frac{\mathrm{d}y }{\mathrm{d} x} = \frac{\mathrm{d} }{\mathrm{d} x}y\\ \frac{\mathrm{d} }{\mathrm{d} x}\frac{\mathrm{d}y }{\mathrm{d} x}+4\frac{\mathrm{d}y }{\mathrm{d} x}+4y=e^x \end{cases}$

此时这是一个关于dy/dx和y的一个微分方程组

令u=(y,dy/dx),我们有

def f(x,u):
	y,dydx = u
	d2ydx2 = np.exp(x)-4*dydx-4*y
	dudx = [dydx,d2ydx2]
	return dudx
u0=[0,0]
x = np.linspace(0,10,100)
sol = odeint(f,u0,x,tfirst=True)

对于更高阶的微分方程也是同样处理

此时我们在回过头来看最初的问题，可以轻松的得到

def f(t,r,k,g):
#k = B*q/m
    x,y,dxdt,dydt = r
    d2xdt2 = -k*dydt
    d2tdt2 = -g + k*dxdt
    return [dxdt,dydt,d2xdt2,d2ydt2]

定义一些常量?

t = np.linspace(0,15,1000)
k = 1
g = 9.8

使用odeint方法

r0 = [0,0,0,0]
sol = odeint(f,r0,t,tfirst = True,args=(k,g))

将图像画出

plt.plot(sol[:,0],sol[:,1])
plt.grid()
plt.show()

odeint解这个方程也十分精确，与sympy相差无几。

最后我们还可以让他实现动画效果

def update_points(num):
    point_ani.set_data(s[:,0][num], s[:,1][num])
    return point_ani,
plt.xlabel('X')
plt.ylabel('Y')
fig,ax = plt.subplots()
plt.plot(s[:,0],s[:,1])
point_ani, = plt.plot(s[:,0][0], s[:,1][0], "ro")
plt.grid(ls="--")
# 生成动画
ani = animation.FuncAnimation(fig, update_points,range(1, len(t)), interval=5, blit=True)
plt.show()