利用Python进行数据分析（解决微分方程）

微分方程的符号解

import numpy
import sympy

y = sympy.symbols("y", cls=sympy.Function)
x = sympy.symbols("x")
eq = sympy.Eq(y(x).diff(x, 2) + 2*y(x).diff(x, 1) + y(x), x**2)
result = sympy.dsolve(eq, y(x))
print(result)

一阶微分方程数值解

'''
dy/dt = sin(t**2)
y(-10) = 1
'''

import scipy
import matplotlib.pyplot as plt

def dy_dt(y, t):    #对于odeint方法默认要求参数是先y后t
    return np.sin(t**2)

y0 = [1]    #y0表示函数的初始值，即t[0]的值
t = np.linspace(-10, 10, 1000)  #函数的取值范围的点
y = scipy.integrate.odeint(dy_dt, y0, t)    #这个函数返回一个在给定的所有点上的y值,形状为 (len(t),len(y0)

plt.plot(t, y)
plt.title("differential equation")

plt.show()

多元一阶微分方程

'''洛伦兹方程组
dx_dt = p*(y-x)
dy_dt = x*(r-z) - y
dz_dt = x*y - b*z
'''

import scipy
import matplotlib.pyplot as plt


def lorenz(w, t, p, r, b):  #w=[x, y, z]
    x, y, z = w
    dx_dt = p*(y-x)
    dy_dt = x*(r-z) - y
    dz_dt = x*y - b*z
    return np.array([dx_dt, dy_dt, dz_dt])

t = np.arange(0, 30, 0.01)  #设置时间
prb = (10, 28, 3)   #设置lorenz里面的参数p，r，b
#求解微分方程
w1 = (0, 1, 0)
slt1 = scipy.integrate.odeint(lorenz, w1, t, args=prb)   #指定args来传递函数参数
w2 = (0, 1.01, 0)
slt2 = scipy.integrate.odeint(lorenz, w2, t, args=prb)  #这里返回的是每个t点对应的[x,y,z]解

fig = plt.figure()

x1 = slt1[:, 0]
y1 = slt1[:, 1]
z1 = np.array([slt1[:, 2]])
x2 = slt2[:, 0]
y2 = slt2[:, 1]
z2 = np.array([slt2[:, 2]])

ax_1 = plt.subplot(projection="3d")
ax_1.plot_wireframe(x1, y1, z1, color='red')
ax_1.set_title("begin with 0,1,0")

ax_2 = plt.subplot(projection="3d")
ax_2.plot_wireframe(x2, y2, z2)
ax_2.set_title("begin with 0,1.01,0")

plt.show()

高阶微分方程

'''
U(t).deff(t, 2) + 2aU(t).deff(t, 1) + (w**2)U(t) = 0
U(0) = U0
U(0).deff(t, 1) = 0
'''
#另V = U(t).diff(t, 1)，可以将方程变形
'''
V = U(t).diff(t, 1)     -->U(t).diff(t, 1) = V
V(t).diff(t, 1) + 2aV(t) + (w**2)U(t) = 0       -->V(t).diff(t, 1) = -2aV(t) - (w**2)U(t)相当于把一个高阶的微分方程编程了多个一阶的微分方程
U(0) = U0
V(0) = 0
'''


import scipy
import matplotlib.pyplot as plt

def deriv(y, t, a, w):
    U, V = y
    du_dt = V
    dv_dt = -2*a*V - (w**2)*U
    return np.array([du_dt, dv_dt])

t = np.arange(0, 20, 0.01)

Y0 = (1, 0)
aw1 = (1, 0.6)  #初始值为1,0.6时（电路上的过阻尼）
slt1 = scipy.integrate.odeint(deriv, Y0, t, args=aw1)
aw2 = (1, 1) #电路中说的临界阻尼
slt2 = scipy.integrate.odeint(deriv, Y0, t, args=aw2)
aw3 = (0.3, 1)   #电路中说的欠阻尼
slt3 = scipy.integrate.odeint(deriv, Y0, t, args=aw3)

plt.plot(t, slt1[:, 0], color='red', label="Over damping")
plt.plot(t, slt2[:, 0], color='green', label="Critical damping")
plt.plot(t, slt3[:, 0], color='blue', label="Underdamping")
plt.legend()

plt.show()