本文介绍如何使用简单的欧拉法求解微分方程,大部分内容出自吴一东老师在他的B站个人空间发布的课程
方法介绍
对于一个一般的微分方程:
{ d y d t = f ( y ( t ) , t ) y ( 0 ) = y 0 \begin{cases} \begin{aligned} \frac{\mathrm{d} y}{\mathrm{d} t} &= f(y(t), t)\\ y(0) &= y_0 \\ \end{aligned} \end{cases} ⎩
⎨
⎧dtdyy(0)=f(y(t),t)=y0
假如我们很难得到他的解析解或者不存在解析解,那我们可以尝试使用欧拉法将微分方程利用泰勒展开构造成一个递推式,通过程序将所有的点求出,进而得到他的数值解。以下说明求法:
-
区间划分
对于一个微分方程我们只关心函数的一段区域内的数值解。
t ∈ [ 0 , T ] t \in [0, T] t∈[0,T] 将区间划成 N N N段,则有 Δ t = T N \Delta t = \frac{T}{N} Δt=NT
t 1 = 0 , t 2 = Δ t , . . . , t n = ( n − 1 ) Δ t , . . . , t N + 1 = T t_1 = 0, t_2 = \Delta t, …, t_{n} = (n-1)\Delta t ,…,t_{N+1} = T t1=0,t2=Δt,…,tn=(n−1)Δt,…,tN+1=T
-
泰勒展开
y ( t + Δ t ) = y ( t ) + d y d t ( t ) Δ t + O ( Δ t 2 ) = y ( t ) + f ( y ( t ) , t ) Δ t + O ( Δ t 2 ) \begin{aligned} y(t+\Delta t) &= y(t) + \frac{\mathrm{d} y} {\mathrm{d} t}(t)\Delta t + O(\Delta t^2) \\ &=y(t) + f(y(t), t)\Delta t + O(\Delta t^2) \end{aligned} y(t+Δt)=y(t)+dtdy(t)Δt+O(Δt2)=y(t)+f(y(t),t)Δt+O(Δt2) -
得到递推公式
y ( t n + 1 ) = y ( t n + Δ t ) = y ( t n ) + f ( y ( t n ) , t n ) Δ t + O ( Δ t 2 ) \begin{aligned} y(t_{n+1}) &= y(t_n + \Delta t)\\ &=y(t_n) + f(y(t_n), t_n)\Delta t + O(\Delta t^2) \end{aligned} y(tn+1)=y(tn+Δt)=y(tn)+f(y(tn),tn)Δt+O(Δt2) -
近似计算
y ( t n ) → t n y n + 1 = y n + f ( y n , t ) Δ t y(t_n) \to t_n\\ y_{n + 1} = y_n +f(y_n, t)\Delta t y(tn)→tnyn+1=yn+f(yn,t)Δt
示例
例1:
{ f ( y , t ) = y d y d t = y y ( 0 ) = 1 \begin{cases} \begin{aligned} f(y, t) &= y\\ \frac{\mathrm{d} y}{\mathrm{d} t} &= y\\ y(0) &= 1 \end{aligned} \end{cases} ⎩
⎨
⎧f(y,t)dtdyy(0)=y=y=1
由题目其实我们很容易得知, y = e t y = e^t y=et 所以我们可以选择解析解对数值解方法进行验证
按上文求解方法我们求其数值解
y n + 1 = y n + y n Δ t y_{n + 1} = y_n + y_n\Delta t yn+1=yn+ynΔt
求解代码如下:
import numpy as np
import matplotlib.pyplot as plt
# 设置初始条件
T = 5
N = 10000
dt = T / N
t = np.linspace(0, T, N + 1)
# 边界条件
y = np.zeros(N + 1)
y[0] = 1
for i in range(0, N):
y[i + 1] = y[i] + y[i] * dt
ex = np.exp(t)
plt.plot(t, y, color='red')
plt.plot(t, ex, color='blue')
plt.show()
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/38225.html