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)。