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