# 微分方程回顾

$f(\frac{d^ny}{dx^n}, \frac{d^{(n-1)}y}{dx^{(n-1)}} ,...,\frac{dy}{dx}, y) =p(x) .$

• 非齐次二级线性微分方程：$\frac{d^2y}{dx^2} +x\frac{dy}{dx}+xy = x+1$
• 谐振子齐次二阶常系数微分方程：
$\frac{d^2u}{dx^2} +\omega ^2u = 0$
• 单摆二阶非线性微分方程：$L\frac{d^2u}{dx^2} +g sin(u) = 0$
• 一阶非线性微分方程：$\frac{dy}{dx}+y^2 = x+1$

# 微分方程：python 解析解（SymPy）

$\frac{d^2x(t)}{dt^2} + 2\gamma \omega _0\frac{dx(t)}{dt} + \omega _0^2 x(t) = 0, \\ initial\ conditions . \begin{cases} x(0) = 0, \\ \frac{dx(t)}{dt}|_{t=0} = 0. \end{cases}$

import numpy as np
from scipy import integrate
#import matplotlib.pyplot as plt
import sympy

def apply_ics(sol, ics, x, known_params):
"""
Apply the initial conditions (ics), given as a dictionary on
the form ics = {y(0): y0, y(x).diff(x).subs(x, 0): yp0, ...},
to the solution of the ODE with independent variable x.
The undetermined integration constants C1, C2, ... are extracted
from the free symbols of the ODE solution, excluding symbols in
the known_params list.
"""

free_params = sol.free_symbols - set(known_params)
## free parameters like C1,C2, that needed to be figured out with ics and boundary conditions
eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics) for n in range(len(ics))]
## form equations with general solution(sol), by substituting the x with 0 and [y(0),dy/dx(x=0)...] with ics.
sol_params = sympy.solve(eqs, free_params)
## solve the algebraic equation set  to get the free parameters expression
## return the solution after substituting the free parameters

return sol.subs(sol_params)

sympy.init_printing()
## initialize the print environment
t, omega0, gamma = sympy.symbols("t, omega_0, gamma", positive=True)
## symbolize the parameters and they can only be positive
x = sympy.Function('x')
## set x to be the differential function, not the variable
ode = x(t).diff(t, 2) + 2 * gamma * omega0 * x(t).diff(t) + omega0**2*x(t)
ode_sol = sympy.dsolve(ode)
## use diff() and dsolve to get the general solution
ics = {x(0): 1, x(t).diff(t).subs(t, 0): 0}
## apply dict to the initial conditions
x_t_sol = apply_ics(ode_sol, ics, t, [omega0, gamma])
## get the solution with ics by calling function apply_ics.
sympy.pprint(x_t_sol)
## pretty print
## solution is as follows
⎛       _______   _______⎞
⎛        γ         1⎞  ω₀⋅t⋅⎝-γ - ╲╱ γ - 1 ⋅╲╱ γ + 1 ⎠   ⎛      γ
x(t) = ⎜- ───────────── + ─⎟⋅ℯ                                + ⎜─────────────
⎜       ________   2⎟                                    ⎜     ________
⎜      ╱  2         ⎟                                    ⎜    ╱  2
⎝  2⋅╲╱  γ  - 1     ⎠                                    ⎝2⋅╲╱  γ  - 1

⎛       _______   _______⎞
1⎞  ω₀⋅t⋅⎝-γ + ╲╱ γ - 1 ⋅╲╱ γ + 1 ⎠
+ ─⎟⋅ℯ
2⎟
⎟
⎠



1. [https://vlight.me/2018/05/01/Numerical-Python-Ordinary-Differential-Equations/]
2. [https://docs.sympy.org/latest/search.html?q=]

# 微分方程：python数值解(SciPY)

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

def plot_direction_field(x, y_x, f_xy, x_lim=(-5, 5), y_lim=(-5, 5), ax=None):
f_np = sympy.lambdify((x, y_x), f_xy, 'numpy')
x_vec = np.linspace(x_lim[0], x_lim[1], 20)
y_vec = np.linspace(y_lim[0], y_lim[1], 20)

if ax is None:
_, ax = plt.subplots(figsize=(4, 4))

dx = x_vec[1] - x_vec[0]
dy = y_vec[1] - y_vec[0]

for m, xx in enumerate(x_vec):
for n, yy in enumerate(y_vec):
Dy = f_np(xx, yy) * dx
Dx = 0.8 * dx**2 / np.sqrt(dx**2 + Dy**2)
Dy = 0.8 * Dy*dy / np.sqrt(dx**2 + Dy**2)
ax.plot([xx - Dx/2, xx + Dx/2], [yy - Dy/2, yy + Dy/2], 'b', lw=0.5)

ax.axis('tight')
ax.set_title(r"$%s$" %(sympy.latex(sympy.Eq(y_x.diff(x), f_xy))), fontsize=18)

return ax

x = sympy.symbols('x')
y = sympy.Function('y')
f = x-y(x)**2

f_np = sympy.lambdify((y(x), x), f)
## put variables (y(x), x) into lambda function f.
y0 = 1
xp = np.linspace(0, 5, 100)
yp = integrate.odeint(f_np, y0, xp)
## solve f_np with initial conditons y0, and x ranges as xp.
xn = np.linspace(0, -5, 100)
yn = integrate.odeint(f_np, y0, xm)

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
plot_direction_field(x, y(x), f, ax=ax)
## plot direction field of function f
ax.plot(xn, yn, 'b', lw=2)
ax.plot(xp, yp, 'r', lw=2)
plt.show()



# 微分方程组：python数值解

$\begin{cases} \frac{dx}{dt} = p(y-x), \\ \frac{dy}{dt} = x(r-z), \\ \frac{dz}{dt} = xy-bz,. \end{cases}$

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()


08-10

01-16 950

06-30 151

10-07 1491

06-13 6688

04-21 2564

10-14 1296

06-19 4808

09-01 4304

05-06 1万+

#### python：使用scipy求解常微分方程

©️2020 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、C币套餐、付费专栏及课程。