微分方程回顾

微分方程是用来描述某一类函数与其导数之间的关系的方程，其解是一个符合方程的函数。（解为常数值的是代数方程）

微分方程的求解是研究微分方程的重点之一，例如解是否存在，存在是否唯一等等。只有少数类型的微分方程存在解析解；无法求得解析解时，可以求数值解（用程序做数值分析），以此确认其解的部分性质。

微分方程按自变量的个数可分为常微分方程（Ordinary Differential Equations）和偏微分方程（Partial Differential Equations），前者只有一个自变量，表达通式形如
f(dnydxn,d(n−1)ydx(n−1),...,dydx,y)=p(x).f(\frac{d^ny}{dx^n}, \frac{d^{(n-1)}y}{dx^{(n-1)}} ,...,\frac{dy}{dx}, y) =p(x) .f(dxndny,dx(n−1)d(n−1)y,...,dxdy,y)=p(x).
后者有两个及以上自变量，例如∂ϕ∂x+x∂ϕ∂y=0.\frac{\partial \phi}{\partial x} + x\frac{\partial \phi}{\partial y} = 0.∂x∂ϕ+x∂y∂ϕ=0.

以常微分方程为例讲一下常见的一些概念，其中包括 （非）齐次，常（变）系数，（非）线性 。

如果第一个公式里左侧为线性微分方程（前提条件），且右侧的p(x)=0p(x)=0p(x)=0，则为齐次线性微分方程，否则非齐次。出现齐次概念的另一个场合是一阶微分方程f(x,y)dy=g(x,y)dxf(x,y) dy = g(x,y)dxf(x,y)dy=g(x,y)dx，具体可以参见Wikipedia 齐次微分方程。简单来说，如果方程的解是齐次函数，那么这个方程就是齐次方程。

常系数和变系数就看函数及其各阶导数前系数是否为常数。

线性，则取决于函数本身是否线性以及函数是否与其导数有乘积，跟自变量无关。
例如：

非齐次二级线性微分方程：d2ydx2+xdydx+xy=x+1\frac{d^2y}{dx^2} +x\frac{dy}{dx}+xy = x+1 dx2d2y+xdxdy+xy=x+1
谐振子齐次二阶常系数微分方程：
d2udx2+ω2u=0\frac{d^2u}{dx^2} +\omega ^2u = 0 dx2d2u+ω2u=0
单摆二阶非线性微分方程：Ld2udx2+gsin(u)=0L\frac{d^2u}{dx^2} +g sin(u) = 0 Ldx2d2u+gsin(u)=0
一阶非线性微分方程：dydx+y2=x+1\frac{dy}{dx}+y^2 = x+1 dxdy+y2=x+1

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

具备解析解的ODE，我们可以利用python的三方包SymPy进行求解。
以求解阻尼谐振子的二阶ODE为例，其表达式为
d2x(t)dt2+2γω0dx(t)dt+ω02x(t)=0,initialconditions.{x(0)=0,dx(t)dt∣t=0=0.\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} dt2d2x(t)+2γω0dtdx(t)+ω02x(t)=0,initial conditions.{x(0)=0,dtdx(t)∣t=0=0.
完整的求解过程如下代码所示。

import numpy as np
from scipy import integrate
#import matplotlib.pyplot as plt
import sympydef apply_ics(sol, ics, x, known_params):"""Apply the initial conditions (ics), given as a dictionary onthe 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 extractedfrom the free symbols of the ODE solution, excluding symbols inthe 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 conditionseqs = [(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 parametersreturn 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⎟                                 ⎟                                 ⎠

参考

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

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

当ODE无法求得解析解时，可以用scipy中的integrate.odeint求数值解来探索其解的部分性质，并辅以可视化，能直观地展现ODE解的函数表达。
以如下一阶非线性（因为函数y幂次为2）ODE为例， dydx=x−y(x)2.\frac{dy}{dx} = x - y(x)^2. dxdy=x−y(x)2.
现用odeint 求其数值解，代码与可视化如下。

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
import sympydef 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) * dxDx = 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 axx = sympy.symbols('x')
y = sympy.Function('y')
f = x-y(x)**2f_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()

右侧曲线跟方向场贴合一致，但左侧蓝线的数值解诡异，明显不满足ODE的表达式，这也说明数值解的局限性。

微分方程组：python数值解

以求解洛伦兹曲线为例，以下方程组代表曲线在xyz三个方向上的速度，给定一个初始点，可以画出相应的洛伦兹曲线。
{dxdt=p(y−x),dydt=x(r−z),dzdt=xy−bz,.\begin{cases} \frac{dx}{dt} = p(y-x), \\ \frac{dy}{dt} = x(r-z), \\ \frac{dz}{dt} = xy-bz,. \end{cases} ⎩⎪⎨⎪⎧dtdx=p(y−x),dtdy=x(r−z),dtdz=xy−bz,.
代码与图如下。

import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as pltdef dmove(Point, t, sets):"""point：present location indexsets：super parameters"""p, r, b = setsx, y, z = Pointreturn 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.01fig = plt.figure()
ax = Axes3D(fig)
ax.plot(P1[:, 0], P1[:, 1], P1[:, 2])
ax.plot(P2[:, 0], P2[:, 1], P2[:, 2])
plt.show()

洛伦兹曲线很好地展示了非线性的直观形态，同时也是展示混沌系统的经典例子。这个例子告诉我们，混沌系统可以是一个确定系统，不一定是随机过程，同时它有初值敏感的特征。

参考
scipy.integrate.odeint用法

Python解微分方程相关推荐

python求解析解,Python解微分方程
Python解微分方程微分方程回顾微分方程:python 解析解(SymPy) 微分方程:python数值解(SciPY) 微分方程组:python数值解微分方程回顾微分方程是用来描述某一类函 ...
python解压zip文件_python-29 python解压压缩包的几种方法
这里讨论使用Python解压例如以下五种压缩文件: .gz .tar .tgz .zip .rar 简单介绍 gz: 即gzip.通常仅仅能压缩一个文件.与tar结合起来就能够实现先打包,再压缩. t ...
delphi dbgrideh 遍历每一个单元格_用Python解数独[1]：求每个单元格的行值域
目录用Python解数独[0] 用Python解数独[1]:求每个单元格的行值域用Python解数独[2]:求列值域和九宫格值域用Python解数独[3]:求总值域用Python解数独[4]: ...
python压缩包怎么安装-详解python解压压缩包的五种方法
这里讨论使用Python解压例如以下五种压缩文件: .gz .tar .tgz .zip .rar 简单介绍 gz: 即gzip.通常仅仅能压缩一个文件.与tar结合起来就能够实现先打包,再压缩. t ...
Matlab符号运算 - 解微分方程
什么是微分方程?如下: 如下图,也是一个微分方程: matlab使用 dsolve 命令解微分方程: 下面来解 dy/dx=4*x^3 这个方程: 一阶,二阶,多阶,都可以解:
python 画图_用python解九宫格以及画图
用python解九宫格的思路很简单,一个是画图部分,用的是turtle库. 演示图像这个九宫格,首先就是画单独的方型,这个函数要自己写: import turtle as t t.speed(0) ...
python deepcopy函数_用Python解数独[6]：递归获得最终答案
目录用Python解数独[0] 用Python解数独[1]:求每个单元格的行值域用Python解数独[2]:求列值域和九宫格值域用Python解数独[3]:求总值域用Python解数独[4]: ...
matlab 向前欧拉公式,向前欧拉公式在Matlab解微分方程初值解的问题
向前欧拉公式在Matlab解微分方程初值解的问题0 fuqilin1202013.07.04浏览527次分享举报用向前欧拉公式(10.8)求解初值问题,dy/dx=-3x+8x-7,y(0)=1,分 ...
matlab中函数或变量无法识别怎么办_用MATLAB巧解微分方程实例分析
点"考研竞赛数学"↑可每天"涨姿势"哦! MATLAB巧解微分方程实例分析王少华西安电子科技大学微分方程求解难, 字母一堆看着烦. 写错数字一时爽, 一直 ...
奥数 python_奥数赛事china夺得冠军！简单思路用Python解经典数学题
2019年第60届国际数学奥林匹克竞赛(IMO)结果出炉,奥数大赛中国夺冠 ,中美两国同时以227的总分并列团体冠军.中国队王者归来!6名队员全部摘金,总成绩荣获世界第一! 不要错过免费学习Pyth ...

Python解微分方程