微分方程求解草稿 未完成 仅为个人学习笔记
Posted 海轰Pro
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了微分方程求解草稿 未完成 仅为个人学习笔记相关的知识，希望对你有一定的参考价值。
目录
微分方程求解
两种方法：解析解和数值解
- 解析解：可以理解为使用数学方式求出微分方程的准确解，求出来就是一串标准的式子
- 数值解：有一些的微分方程解析解不好求或者求不出，这时候可以求它的数值解，就是使用一些方法模拟出标准解的一个大致图像（逼近解析解）
以上为个人理解 如有错误 请您指出～
例 - 1
比如求下列微分方程的解：
{ d y d x = 2 x + 2 y ( 0 ) = 1 \\begin{cases} \\frac {dy}{dx} = 2x + 2\\\\ y(0) = 1 \\end{cases} {dxdy=2x+2y(0)=1
标准解
其实可以非常轻松的求出答案： f ( x ) = x 2 + 2 x + 1 f(x) = x^2 + 2x + 1 f(x)=x2+2x+1
f ‘ ( x ) = 2 x + 2 f^`(x) = 2x + 2 f‘(x)=2x+2 这个就不解释了哈
绘制 f ( x ) = x 2 + 2 x + 1 f(x) = x^2 + 2x + 1 f(x)=x2+2x+1 的图像
这个就是此题微分方程的一个标准解（解析解），观察其图像
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,5,100)
y = x**2 + 2*x + 1
plt.plot(x,y)
plt.show()
求解析解
使用sympy求解析解
也就是直接求出 f ( x ) = x 2 + 2 x + 1 f(x) = x^2 + 2x + 1 f(x)=x2+2x+1 并显示
首先我们对原来的式子进行一个简单的变形
{ d y d x − 2 x − 2 = 0 ( 保 证 等 式 右 边 为 0 ) y ( 0 ) = 1 \\begin{cases} \\frac {dy}{dx} - 2x - 2 = 0 (保证等式右边为0)\\\\ y(0) = 1 \\end{cases} {dxdy−2x−2=0(保证等式右边为0)y(0)=1
也可以写为
{ f ‘ ( x ) − 2 x − 2 = 0 ( 保 证 等 式 右 边 为 0 ) y ( 0 ) = 1 \\begin{cases} f^`(x) - 2x - 2 = 0 (保证等式右边为0)\\\\ y(0) = 1 \\end{cases} {f‘(x)−2x−2=0(保证等式右边为0)y(0)=1
import sympy
def apply_ics(sol, ics, x, known_params):
free_params = sol.free_symbols - set(known_params)
eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics) for n in range(len(ics))]
sol_params = sympy.solve(eqs, free_params)
return sol.subs(sol_params)
# 初始化打印环境
sympy.init_printing()
# 标记参数，且均为正
x = sympy.symbols("x", positive = True)
# 标记y是微分函数，非变量
y = sympy.Function("y")
# ode 微分方程等号左边的部分，等号右边为0
ode = y(x).diff(x) - 2*x - 2
# 求解
ode_sol = sympy.dsolve(ode)
# 初始条件：字典匹配
ics = {y(0): 1}
# 显示
x_t_sol = apply_ics(ode_sol, ics, x, [])
sympy.pprint(x_t_sol)
运行结果
可以发现，求出了该微分方程的一个非常标准的解
求数值解
使用integrate.odeint求数值解
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
# 引入中文字体
font = {
"family": "Microsoft YaHei"
}
matplotlib.rc("font", **font)
def dy_dt(y, x): # 定义导数函数 f(y,x)
return 2*x + 2
y0 = [1] #
x = np.arange(0,5,0.01) # (start,stop,step) start = t0
y = odeint(dy_dt, y0, x) # 求解微分方程初值问题
# 绘图
plt.plot(x, y)
plt.title("f(x)的数值解(近似解)")
plt.show()
运行结果
可以发现数值解得出的结果与标准解图像非常逼近，在一定程度上可以作为标准解来使用、观察
Notes
简单理解，解析解就是使用纯数学方法来求微分方程的解，得到的结果是一个非常标准的微分方程解；而数值解则是在某一个区间内对解的图像形状进行一种模拟、逼近，得到的是在一个区间内标准解的一个近似图像
例 - 2
求下列微分方程的解
{
(
1
+
x
2
)
y
′
+
2
x
y
=
x
e
x
2
y
(
0
)
=
−
1
2
\\begin{cases} (1+x^2)y^{'} + 2xy = xe^{x^2}\\\\ y(0) = - \\frac{1}{2} \\end{cases}
{(1+x2)y′+2xy=xex2y(0)=−21
标准解
运用数学方法可以求出标准解为：
y
=
e
x
2
−
2
2
(
1
+
x
2
)
y = \\frac {e^{x^2 } -2 }{2(1+x^2)}
y=2(1+x2)ex2−2
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,3,100)
y = (np.e**(x**2) - 2) / (2*(1+x**2))
plt.plot(x,y)
plt.show()
求解析解
改写一下
{
(
1
+
x
2
)
y
′
+
2
x
y
−
x
e
x
2
=
0
（
等
式
右
边
必
须
为
0
）
y
(
0
)
=
−
1
2
\\begin{cases} (1+x^2)y^{'} + 2xy - xe^{x^2} = 0（等式右边必须为0）\\\\ y(0) = - \\frac{1}{2} \\end{cases}
{(1+x2)y′+2xy−xex2=0（等式右边必须为0）y(0)=−21
import sympy
def apply_ics(sol, ics, x, known_params):
free_params = sol.free_symbols - set(known_params)
eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics) for n in range(len(ics))]
sol_params = sympy.solve(eqs, free_params)
return sol.subs(sol_params)
# 初始化打印环境
sympy.init_printing()
# 标记参数，且均为正
x = sympy.symbols("x", positive = True)
# 标记y是微分函数，非变量
y = sympy.Function("y")
# ode 微分方程等号左边的部分，等号右边为0
ode = (1+x**2)*y(x).diff(x,1) + 2*x*y(x) - x*np.e**(x**2)
# 求解
ode_sol = sympy.dsolve(ode)
# 初始条件：字典匹配
ics = {y(0): -0.5}
# 显示
x_t_sol = apply_ics(ode_sol, ics, x, [])
sympy.pprint(x_t_sol)
与标准解完全一致（分子分母同时乘以2就行）
求数值解
改写一下
{
y
′
=
x
e
x
2
−
2
x
y
1
+
x
2
y
(
0
)
=
−
1
2
\\begin{cases} y^{'}= \\frac {xe^{x^2} - 2xy}{1+x^2}\\\\ y(0) = - \\frac{1}{2} \\end{cases}
{