Реализация метода Эйлера в Python дает стабильный результат, но он должен быть нестабильным

Я пытаюсь решить это дифференциальное уравнение с помощью метода Эйлера с использованием 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):

введите здесь описание изображения


comment
Не используйте math.e**(-10*x), используйте np.exp(-10*x), так как a^b все равно реализовано как exp(b*log(a)).   -  person Lutz Lehmann    schedule 08.05.2020


Ответы (2)


В зависимости от размера временного шага интегрирование Эйлера может быть устойчивым или неустойчивым, поскольку это явный метод. Вы выбрали довольно маленький временной шаг. Если вы увеличите его, вы начнете видеть колебания, как показано на рисунке ниже.

Колебания

Вот небольшая тестовая программа, которую я написал (попробуйте медленно увеличивать переменную 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()
person Sourabh Bhat    schedule 28.09.2018
comment
Привет, спасибо за попытку помочь мне! Я сделал то, что вы сказали, и увеличил свой шаг до 0,9. Я также увеличил интервал от 0 до 50, чтобы лучше видеть. Однако моя реализация не генерирует колебания! Он генерирует монотонную ошибку вниз. У вас есть понимание, почему это происходит? Спасибо за ваш код. Но мне бы очень хотелось понять, почему моя реализация терпит неудачу. - person Pedro Delfino; 28.09.2018
comment
В вашем коде также для delta_x = 0,22 и выше я вижу видимые колебания. Я использую Python 3. - person Sourabh Bhat; 30.09.2018
comment
Ах, да! Ты прав. Колебания на самом деле *происходят в моей реализации, просто они были довольно маленькими, я их не видел! Спасибо чувак. - person Pedro Delfino; 30.09.2018

Это классическая задача численного анализа. Метод Эйлера нестабилен только в том случае, если размер вашего шага слишком велик. Точное решение вашего дифференциального уравнения:

у(т) = ехр(-10т)*у(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.

Неявный метод может позволить вам обойти это ограничение временного шага.

person numerics dude    schedule 08.05.2020