Python 中的 Euler 方法实现给出了稳定的结果，但它应该是不稳定的

Question

我正在尝试使用欧拉方法求解这个微分方程 Python3:

根据 Wolfram Alpha，这是正确方程的绘图。

同样，根据 Wolfram Alpha，在这种情况下，经典的欧拉方法应该不会稳定，正如您在区间结束时看到的那样：

然而，在我的实现中，欧拉方法提供了一个稳定的结果，很奇怪。我想知道我的实施由于某种原因是错误的。尽管如此，我还是找不到错误。

我生成了一些点和绘图，比较了我的近似值和函数的分析输出。蓝色为对照组的分析结果。红色部分是我的实施输出：

这是我的代码：

import math
import numpy as np
from matplotlib import pyplot as plt
import pylab

def f(x):

    return (math.e)**(-10*x)

def euler(x):

    y_init = 1
    x_init = 0

    old_dy_dx = -10*y_init

    old_y = y_init 

    new_y = None

    new_dy_dx = None

    delta_x = 0.001

    limite = 0

    while x>limite:

        #for i in range(1,6):

        new_y = delta_x*old_dy_dx + old_y
        #print ("new_y", new_y)

        new_dy_dx = -10*new_y
        #print ("new dy_dx", new_dy_dx)

        old_y = new_y
        #print ("old_y", old_y)

        old_dy_dx = new_dy_dx
        #print ("old delta y_delta x", old_dy_dx)
        #print ("iterada",i)

        limite = limite +delta_x

    return new_y

t = np.linspace(-1,5, 80)

lista_outputs = []

for i in t:
    lista_outputs.append(euler(i))
    print (i)

# red dashes, blue squares and green triangles
plt.plot(t, f(t), 'b-', label='Output resultado analítico')
plt.plot(t , lista_outputs, 'ro', label="Output resultado numérico")
plt.title('Comparação Euler/Analítico - tolerância: 0.3')
pylab.legend(loc='upper left')
plt.show()

感谢您的帮助。

============================================= ===============

更新

在@SourabhBhat 的帮助下，我能够看到我的实施实际上是正确的。确实，这造成了不稳定。除了增加步长外，我还需要放大才能看到它的发生。

下图说明了一切（步长为 0.22）：

Answer 1

根据时间步长，欧拉积分可以是稳定的也可以是不稳定的，因为它是一种显式方法。您选择了一个非常小的时间步长。如果增加它，您将开始看到振荡，如下图所示。

这是我写的一个小测试程序（尝试慢慢增加 steps 变量 [20,30,40,50....]）：

import numpy as np
import matplotlib.pyplot as plt

steps = 20


def exact_solution(t):
    return np.exp(-10.0 * t)


def numerical_solution(y0, dt, num_steps):
    y = np.zeros(num_steps + 1)
    y[0] = y0
    for step in range(num_steps):
        y[step + 1] = y[step] - 10.0 * y[step] * dt

    return y


if __name__ == "__main__":
    t0 = 0
    time = np.linspace(t0, 5, steps + 1)
    num_sol = numerical_solution(exact_solution(t0), time[1] - time[0], steps)
    exact_sol = exact_solution(time)

    plt.plot(time, num_sol, ".--", label="numerical")
    plt.plot(time, exact_sol, label="exact")
    plt.legend(loc="best")
    plt.show()

Answer 2

这是一个经典的数值分析问题。仅当您的步长太大时，欧拉方法才不稳定。你的微分方程的精确解是：

y(t) = exp(-10t)*y(0)

那是你的初始数据简单地乘以 exp(-10t)，它总是一个满足 0 < exp(-10t) < 1 的数字。因此你的结果变小了。

单时间步的欧拉方法得出：

y(t + dt) = (1-10dt)*y(t)

因此，您的数据只是乘以 (1-10dt)。您希望 0 < (1-10dt) < 1 来模拟原始 ODE。即：

0 < dt < 1/10

如果你在这个范围内选择你的步长，欧拉法将是稳定的。

较大的 dt 将导致因子 (1-10*dt) 变为负数，并且您的解决方案将在每一步后改变其符号。选择更大的 dt，dt > 2/10 将导致您的解决方案呈指数增长，因为那时 |1-10*dt| > 1.

隐式方法可以让您绕过这个时间步长限制。

Python 中的 Euler 方法实现给出了稳定的结果，但它应该是不稳定的

Euler Method implementation in Python gives a stable result but it should be unstable

numerical-analysis

numerical-methods

python-3.x

differential-equations

numerical-stability