欧拉方法（显式和隐式）

Posted 2021-03-28

我想为以下模型实现Euler的方法（显式和隐式）（https://en.wikipedia.org/wiki/Euler_method）：

x(t)' = q(x_M -x(t))x(t)
x(0) = x_0 

其中q，x_M和x_0是实数。

我已经知道该方法的（理论上）实现。但我无法弄清楚我可以插入/更改模型的位置。有人可以帮忙吗？

编辑：你是对的。我没理解这个方法。现在，几个小时后，我想我真的明白了！使用显式方法，我很确定（尽管如此：任何人都可以查看我的代码吗？）

有了隐式实现，我不太确定它是否正确。可以请任何人看看隐式方法的实现，并给我一个反馈什么是正确/不好？

def explizit_euler():
    ''' x(t)' = q(xM -x(t))x(t)
    x(0) = x0'''
    q = 2. 
    xM = 2
    x0 = 0.5 
    T = 5 
    dt = 0.01 
    N = T / dt 
    x = x0
    t = 0.
    for i in range (0 , int(N)):
        t = t + dt
        x = x + dt * (q * (xM - x) * x)
        print '%6.3f %6.3f' % (t, x)



def implizit_euler():
    ''' x(t)' = q(xM -x(t))x(t)
    x(0) = x0'''
    q = 2.
    xM = 2
    x0 = 0.5
    T = 5
    dt = 0.01
    N = T / dt
    x = x0
    t = 0.
    for i in range (0 , int(N)):
        t = t + dt
        x = (1.0 / (1.0 - q *(xM + x) * x))
        print '%6.3f %6.3f' % (t, x)

答案

优先注意事项：尽管总体思路应该是正确的，但我在编辑框中完成了所有代数，因此可能存在错误。在使用任何非常重要的东西之前，请自行检查。

我不确定你是如何看待“隐含的”公式的

    x = (1.0 / (1.0 - q *(xM + x) * x))

但这是错误的，你可以通过比较你的“明确”和“隐含”结果来检查它们：它们应该略有不同但是使用这个公式它们会大大分歧。

要理解隐式Euler方法，您应该首先了解显式方法背后的想法。这个想法非常简单，并且在维基的Derivation部分进行了解释：由于衍生的y'(x)是(y(x+h) - y(x))/h的极限，你可以将y(x+h)近似为y(x) + h*y'(x)的小h，假设我们的原始微分方程是

y'(x) = F(x, y(x))

请注意，这只是一个近似值而不是精确值的原因是，即使在小范围的[x, x+h]上，导数y'(x)也会略有变化。这意味着如果你想得到更好的近似y(x+h)，你需要在y'(x)范围内更好地逼近“平均”导数[x, x+h]。让我们称之为近似y'。这种改进的一个想法是同时找到y'和y(x+h)说我们想找到y'和y(x+h)，y'实际上是y'(x+h)（即最后的衍生物）。这导致以下方程组：

y'(x+h) = F(x+h, y(x+h))
y(x+h) = y(x) + h*y'(x+h)

这相当于一个“隐含”的等式：

y(x+h) - y(x) = h * F(x+h, y(x+h))

它被称为“隐含”，因为在这里目标y(x+h)也是F的一部分。请注意，在wiki文章的Modifications and extensions部分中提到了非常类似的等式。

所以现在谈谈你的方程式

x(t+dt) - x(t) = dt*q*(xM -x(t+dt))*x(t+dt)

或者等价的

dt*q*x(t+dt)^2 + (1 - dt*q)*x(t+dt) - x(t) = 0

这是一个二次方程，有两个解决方案：

x(t+dt) = [(dt*q - 1) ± sqrt((dt*q - 1)^2 + 4*dt*q*x(t))]/(2*dt*q)

显然，我们希望解决方案“接近”x(t)，即+解决方案。所以代码应该是这样的：

    b = (q * xM * dt - 1)
    x = (b + (b ** 2 + 4 * q * x * dt) ** 0.5) / 2 / q / dt