python解常微分方程（组）

最新推荐文章于 2024-06-04 08:30:00 发布

whr804

最新推荐文章于 2024-06-04 08:30:00 发布

阅读量9.6k

点赞数 27

文章标签： python

本文链接：https://blog.csdn.net/weixin_61268973/article/details/122462611

版权

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

一周前看到了一个视频，于是我便想用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}$

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

#导入
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)

关注