欧拉法（Euler‘s method）求常微分方程(ODE)近似解

最新推荐文章于 2025-03-23 22:07:11 发布

生信小兔

最新推荐文章于 2025-03-23 22:07:11 发布

阅读量1.4w

点赞数 12

分类专栏：数学基础文章标签：算法 python

本文链接：https://blog.csdn.net/weixin_60737527/article/details/125013910

版权

数学基础专栏收录该内容

14 篇文章

订阅专栏

求微分方程的近似解在实际问题的解决中具有极其重要的作用，因为在大多情况下我们很难求出微分方程的函数表达式，或是表达式过于复杂，此时便可以用数值计算的方法求该微分方程的近似解。欧拉法便是其中一种最为简单的求常微分方程近似解的方法，该方法较为简单，同时也有比较明显的缺陷，但也是其它的求解微分方程近似解方法的基础。

算法原理

欧拉法求解过程是一个递归的过程，这个思想和牛顿法、梯度下降法是相似的。并且它将函数离散化，分割成一个个小段来求解。欧拉法求解的常微分方程的形式通常为：

$\begin{Bmatrix}\frac{dy}{dx}=f(x,y) ,&a\leq x\leq b\\y(a)=y_{0} \end{matrix}$

右边是关于x,y的任意函数，如2x,3y+2或 $x^{2}+y^{2}$ 等都是可以的,但y是x的函数，并且函数定义在区间[a,b]内，初值条件为x=a时函数值为 $y_{0}$ 。欧拉法求解思想是把区间[a,b]分成n份，那么每一小段长度就是 $h=\frac{b-a}{n}$ ,初始时左端点a=x0,那么下一个节点x1=x0+h,再下一个就是x2=x1+h，......以此类推，最后 $X_{n}$ =b。我们对上面的微分方程两边同时从 $x_{k}$ 到 $x_{k+1}(k=0,1...n-1)$ 进行积分，就得到下面结果，即：

$y(x_{k+1})-y(x_{k})=\int_{x_{k}}^{x_{k+1}}f(x,y(x))dx$

这样便可以由初始条件y(x0)= $y_{0}$ 与该条件式逐个递推，依次求得y(x1),y(x2)...y(xn)各个函数值。对于右端的积分，在实际计算过程中也需要用近似表达式来替代。欧拉法使用的是由矩形面积来近似替代，即：

$\int_{x_{k}}^{x_{x+1}}f(x,y(x))dx\approx f(x_{k},y(x_{k}))(x_{k+1}-x_{k}) \\ =f(x_{k},y(x_{k}))h \ (k=0,1...n-1)$

这样我们就得到了最终的递归表达式：

$y_{k+1}=y_{k}+f(x_{k},y_{k})h$

从而能够比较简便地进行计算。

算法实现

兔兔这里以一个微分方程为例：

$\begin{Bmatrix}\frac{dy}{dx} =y-\frac{2x}{y},&0\leq x\leq 1 \\y(0)=1\end{}$

对应前面的表达式，这里 $f(x,y(x))=y-\frac{2}{y}$ ,所以计算过程为 $y_{1}=y_{0}+h(\frac{2x_{0}}{y_{0}})$ ,并且x0=0,y0=1，之后再求解y2,y3...yn。h根据需要可以做不同调整，如果h=0.1就是分10份，h越小，分的份数越多，点越来越密集，就可以体现函数的形式。

import numpy as np
import matplotlib.pyplot as plt
y0=1 #初始条件
def f(x,y):
    return y-2*x/y
def Eular(y0,n,f):
    y=[y0] #储存各个点的函数值
    y0=y0
    for x in np.linspace(0,1,n):
        '''求函数在区间[0,1]近似解，分n份'''
        y1=y0+f(x,y0)/n
        y.append(y1)
        y0=y1
    return y
y=Eular(1,10,f=f)
x=np.linspace(0,1,len(y))
plt.plot(x,y,color='blue')
plt.show()

该函数是=的精确解为 $y=\sqrt{1+2x}$ ,我们最终把精确解，n=10,n=100的结果进行比较，图像如下所示：