Я пытаюсь решить это дифференциальное уравнение с помощью метода Эйлера с использованием 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):
math.e**(-10*x)
, используйтеnp.exp(-10*x)
, так какa^b
все равно реализовано какexp(b*log(a))
. - person Lutz Lehmann   schedule 08.05.2020