解常微分方程

解常微分方程问题

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

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

 

sympy中的dsolve方法

Python例程

#导入
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])
print("微分方程:",sol)   #现在打印出sol

###这个式子非常冗杂，用trigsimp()方法化简
x = trigsimp(sol[0].rhs)
y = trigsimp(sol[1].rhs)
print("化简:",x,y)

###将里面的积分常数算出
#定义积分变量，避免报错,观察上式输出式子中有几个C这里就定义几个
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])
print("积分常数：",l)

###将积分常量替换到x和y里面,我们就得到了解的最终形式
x = x.subs(l)
y = y.subs(l)
print("最终形式：",x,y)

#作图
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()

解一阶常微分方程：dy/dx=y
Python例程

from sympy import *
f = symbols('f', cls=Function)#定义函数标识符
x = symbols('x')#定义变量
eq = Eq(diff(f(x),x,1),f(x))#构造等式，即dy/dx=y
#diff(函数,自变量,求导次数)
print(dsolve(eq, f(x)))
pprint(dsolve(eq, f(x)))#以"pretty"形式打印方程的解

解二阶常微分方程：y"+py'+qy=0

Python例程

from sympy import *
f = symbols('f', cls=Function)#定义函数标识符
x,p,q = symbols('x,p,q')#批量定义变量
eq = Eq(diff(f(x),x,2)+p*diff(f(x),x,1)+q*f(x),0)
#构造方程，即y"+py'+qy=0
#diff(函数,自变量,求导次数)
print(dsolve(eq, f(x)))
pprint(dsolve(eq, f(x)))#以"pretty"形式打印方程的解

 

常微分方程绘图

Python例程

#绘图
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def diff(y, x):
   return np.array(1/x)
   # 上面定义的函数在odeint里面体现的就是dy/dx =1/x    #定义常微分方程
x = np.linspace(1, 10, 100)  # 给出x范围，(起始，终值，分割段数)
y = odeint(diff, 0, x)  # 设函数初值为0，即f(1)=0
plt.xlabel('x')
plt.ylabel('y')
plt.title("y=ln(x)")    #定义图的标题
plt.grid()#绘制网格
plt.plot(x, y)  # 将x,y两个数组进行绘图
plt.show()

 

单摆运动：y"+gsin(y)/l=0,g为重力加速度，l为摆长

Python例程

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np

g = 9.8
l = 1
#重力加速度为9.8m/s,摆长为1m，y"+gsin(y)/l=0,g为重力加速度，l为摆长
def diff(d_list, t):#我们可以将一个二阶常微分方程分解为含有两个方程的一阶常微分方程组
   omega, theta = d_list
   return np.array([-g/l*np.sin(theta), omega])
'''
深入剖析diff函数：diff的左边代表dω/dt和dθ/dt，由于函数返回的是数组类型，我们
可以用这种方式构建一个微分方程组：dθ/dt=ω,dω/dt=-gsin(θ)/l
'''
t = np.linspace(0, 10, 1000)
result = odeint(diff, [0, 30/180*np.pi], t)
# odeint中第二个是初始单摆角度30度，无法采用小角近似
plt.xlabel('x')
plt.ylabel('y')
plt.title("dθ/dt=ω,dω/dt=-gsin(θ)/l")
plt.plot(t, result)  # 输出ω和θ随时变化曲线，即方程解
plt.show()

 

洛伦兹曲线：

　　以下方程组代表曲线在xyz三个方向上的速度，给定一个初始点，可以画出相应的洛伦兹曲线。

Python例程

import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def dmove(Point, t, sets):
    """
    point：present location index
    sets：super parameters
    """
    p, r, b = sets
    x, y, z = Point
    return np.array([p * (y - x), x * (r - z), x * y - b * z])

t = np.arange(0, 30, 0.001)
P1 = odeint(dmove, (0., 1., 0.), t, args=([10., 28., 3.],))  #
## (0.,1.,0.) is the initial point; args is the set for super parameters
P2 = odeint(dmove, (0., 1.01, 0.), t, args=([10., 28., 3.],))
## slightly change the initial point from y = 1.0 to be y = 1.01

fig = plt.figure()
ax = Axes3D(fig)
ax.plot(P1[:, 0], P1[:, 1], P1[:, 2])
ax.plot(P2[:, 0], P2[:, 1], P2[:, 2])
plt.show()