python 微分方程组 python求微分方程的解

sympy、numpy、scipy、matplotlib是强大的处理数学问题的库，可以执行积分、求解常微分方程、绘图等功能，其开源免费的优势可以与MATLAB媲美。

• 一阶常微分方程

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"形式打印方程的解

输出结果：

大家可以看出两种打印方式的区别。

• 二阶常微分方程

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"形式打印方程的解

输出结果：

• 对可解的微分方程进行绘图：
  （特别注意odeint的用法）

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()#打印图表

输出结果：

• 不可求出解析解的常微分方程
  我们考虑单摆的运动：θ"+gsinθ/l=0

from sympy import *
f = symbols('f', cls=Function)#定义函数标识符
x,g,l= symbols('x,g,l')#批量定义变量
eq = Eq(diff(f(x),x,2)+g*sin(f(x))/l,0)
#构造单摆运动的标准常微分方程，即y"+gsin(y)/l=0,g为重力加速度，l为摆长
#diff(函数,自变量,求导次数)
print(dsolve(eq, f(x)))
pprint(dsolve(eq, f(x)))#以"pretty"形式打印方程的解

输出结果：

可以看出，dsolve功能对此类常微分方程无能为力。但是，我们可以通过scipy和matplotlib，直观地通过图像呈现其数值解：

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

g = 9.8
l = 1
#重力加速度为9.8m/s,摆长为1m
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()

输出结果：