解常微分方程问题

目前创新互联已为数千家的企业提供了网站建设、域名、网页空间、网站托管运营、企业网站设计、南岳网站维护等服务，公司将坚持客户导向、应用为本的策略，正道将秉承"和谐、参与、激情"的文化，与客户和合作伙伴齐心协力一起成长，共同发展。

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

则可直接列出粒子的运动方程，将这个方程分解成x和y两个方向，联立即可求得该方程组的解。

sympy中的dsolve方法

Python例程

 1 #导入
 2 from sympy import *
 3 import numpy as np
 4 import matplotlib.pyplot as plt
 5 init_printing()
 6 
 7 ###首先声明符号x,y,q,m,B,g
 8 #参量
 9 q,m,B,t,g = symbols('q m B t g',real=True,nonzero=True,positive=True)
10 #函数
11 x,y = symbols('x y',real=True,nonzero=True,positive=True,cls=Function)
12 
13 ###再将微分方程表示出来
14 eq1 = Eq(x(t).diff(t,t),-q*B*y(t).diff(t)/m)
15 eq2 = Eq(y(t).diff(t,t),-g+q*B*x(t).diff(t)/m)
16 sol = dsolve([eq1,eq2])
17 print("微分方程:",sol)   #现在打印出sol
18 
19 ###这个式子非常冗杂，用trigsimp()方法化简
20 x = trigsimp(sol[0].rhs)
21 y = trigsimp(sol[1].rhs)
22 print("化简:",x,y)
23 
24 ###将里面的积分常数算出
25 #定义积分变量，避免报错,观察上式输出式子中有几个C这里就定义几个
26 var('C1 C2 C3 C4')
27 #x(0)=0
28 e1 = Eq(x.subs({t:0}),0)
29 #x'(0)=0,要将subs放在diff后面
30 e2 = Eq(x.diff(t).subs({t:0}),0)
31 #y(0)=0
32 e3 = Eq(y.subs({t:0}),0)
33 #y'(0)=0
34 e4 = Eq(y.diff(t).subs({t:0}),0)
35 l = solve([e1,e2,e3,e4],[C1,C2,C3,C4])
36 print("积分常数：",l)
37 
38 ###将积分常量替换到x和y里面,我们就得到了解的最终形式
39 x = x.subs(l)
40 y = y.subs(l)
41 print("最终形式：",x,y)
42 
43 #作图
44 ts = np.linspace(0,15,1000)
45 consts = {q:1,B:1,g:9.8,m:1}
46 fx = lambdify(t,x.subs(consts),'numpy')
47 fy = lambdify(t,y.subs(consts),'numpy')
48 plt.plot(fx(ts),fy(ts))
49 plt.grid()
50 plt.show()