2

我正在尝试使用 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）：

在此处输入图像描述

4

2 回答 2

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()
于 2018-09-28T05:07:02.973 回答
0

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

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。

隐式方法可以让您绕过此时间步长约束。

于 2020-05-08T09:22:13.437 回答