2023年6月29日发(作者：)

python解常微分⽅程（组）本⼈⽬前初三，能⼒所限，如有不⾜之处，还望多多指教。⼀周前看到了⼀个，于是我便想⽤python来求解这个问题。〇、分析假设在平⾯内有⼀带电粒⼦q，质量为m。空间内存在匀强磁场B，⽅向垂直于平⾯向内即沿z轴负半轴，以及⼀个沿y轴负半轴的重⼒场。带电粒⼦从磁场内O点释放。则可直接列出粒⼦的运动⽅程将这个⽅程分解成x和y两个⽅向联⽴即可求得该⽅程组的解。⼀、sympy中的dsolve⽅法#导⼊from sympy import *import numpy as npimport as pltinit_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)=0e1 = Eq(({t:0}),0)#x'(0)=0,要将subs放在diff后⾯e2 = Eq((t).subs({t:0}),0)#y(0)=0e3 = Eq(({t:0}),0)#y'(0)=0e4 = Eq((t).subs({t:0}),0)l = solve([e1,e2,e3,e4],[C1,C2,C3,C4])l紧接着，将积分常量替换到x和y⾥⾯,我们就得到了解的最终形式x = (l)y = (l)x,y当然，这个解还可以写成这种形式⽤plt画图来看看ts = ce(0,15,1000)consts = {q:1,B:1,g:9.8,m:1}fx = lambdify(t,(consts),'numpy')fy = lambdify(t,(consts),'numpy')(fx(ts),fy(ts))()()但是sympy有⼀个缺点，当微分⽅程很复杂时，它会直接罢⼯。于是另⼀个新的⽅法替我们解决了这个问题⼆、ate中的odeint⽅法#导⼊import numpy as npfrom ate import odeintimport as pltimport ion as animationfrom math import e

我们⾸先要创建⼀个函数，使它可以表⽰这个微分⽅程。对于这个函数应该具备怎样的形式，先从⼀阶微分⽅程开始⽐如说，我们要求解下⾯这个⽅程它的解通过⼀些简单的分离变量即得下⾯⽤odeint来解出它，令y(1)=1#⼀阶def f(x,y): #将⼀阶导数⽤y和x构成的函数来表⽰ dydx = -y/x return dydx#初始条件y0 = 1#初值条件是在⾃变量x的下限取到。如y(1)就要⽣成⼀个从1开始的范围x = ce(1,5,100)(0,5)(0,5)#odeint(）⽅法,tfirst属性指的是函数f的第⼀个参数是⾃变量sol = odeint(f,y0,x,tfirst=True)(x,sol[:,0])()()同理，对于下⾯的⼀阶微分⽅程组此时我们设向量u=(x,y) （列向量），⽅程组化为其实对于任意的⼀阶微分⽅程，都可以写成这样的形式f的含义就是此时微分⽅程就可以化为⽽f(t,u)正是我们要找的那个函数odeint中⽀持向量输⼊,因此，可以这样构造这个函数def f(t,u): # = 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 = ce(0,10,100)sol = odeint(f,u0,t,tfirst=True)因此，对于⼆阶常微分⽅程我们⾸先把这个⽅程化解成这样的形式此时这是⼀个关于dy/dx和y的⼀个微分⽅程组令u=(y,dy/dx),我们有def f(x,u): y,dydx = u d2ydx2 = (x)-4*dydx-4*y dudx = [dydx,d2ydx2] return dudxu0=[0,0]x = ce(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 = ce(0,15,1000)k = 1g = 9.8使⽤odeint⽅法r0 = [0,0,0,0]sol = odeint(f,r0,t,tfirst = True,args=(k,g))将图像画出(sol[:,0],sol[:,1])()()odeint解这个⽅程也⼗分精确，与sympy相差⽆⼏。最后我们还可以让他实现动画效果def update_points(num): point__data(s[:,0][num], s[:,1][num]) return point_ani,('X')('Y')fig,ax = ts()(s[:,0],s[:,1])point_ani, = (s[:,0][0], s[:,1][0], "ro")(ls="--")# ⽣成动画ani = imation(fig, update_points,range(1, len(t)), interval=5, blit=True)()如果要深度学习这两个库的话，还是推荐官⽅⽂档最近还准备在更新⼏篇（如果期末不挂的话） 注：我在b站上也发布了这篇⽂章