不同步长的欧拉方法。如何更改算法的代码以考虑不同的步长值?
Posted
技术标签:
【中文标题】不同步长的欧拉方法。如何更改算法的代码以考虑不同的步长值?【英文标题】:Euler's method for different step sizes . How does one change the code of the algorithm to account for different values of the step size? 【发布时间】:2020-08-08 19:01:38 【问题描述】:我有一个用于数值微分方程问题的算法,称为欧拉法。本质上,欧拉方法近似于微分方程的解。我的函数适用于单步长(值h
),但我正在尝试更改代码以允许我循环3个不同的值h(通过将h
从单个值更改为可能值的列表) .但是,我编写的函数没有充分循环我的值。我是 python 新手,以前使用过 R。有人可以告诉我如何正确地做到这一点。
我的代码适用于步长 h 的单个值是:
from math import exp # exponential function
dy = lambda x,y: x*y
f = lambda x: exp(x**2/2) # analytical solution function
x = 0 # Intial value X_0
xn = 2 # Final Value
y = 1 # value of y(x0)
h = 0.2 # stepsize
n = int((xn-x)/h)
print ('x \t\t y (Euler h=) \t y (analytical)'.format(h))
print ('%f \t %f \t %f'% (x,y,f(x)))
for i in range(n):
y += dy(x, y)*h
x += h
print ('%f \t %f \t %f'% (x,y,f(x)))
x y (Euler h=0.5) y (analytical)
0.000000 1.000000 1.000000
0.500000 1.000000 1.133148
1.000000 1.250000 1.648721
1.500000 1.875000 3.080217
2.000000 3.281250 7.389056
我想将 h 更改为 h=[0.01,0.2,0.5]
并获得值,然后创建显示解析解和欧拉方法解在不同步长值下的图。
如果这是一个简单的问题,我再次道歉。我是 python 编程的新手,经常犯一些错误,下面是我迄今为止最好的尝试。我还没有将我的 x 值存储到容器中,因为我的函数没有循环遍历 h 值。我正在尝试编写一个嵌套的 for 循环,其中外部循环遍历 h 值并存储这些值并将它们绘制为一条线,然后迭代到 h 的第二个值并执行相同的操作,最后这些值可以是放在一个地块上。
# Improved to allow plotting different values
import matplotlib.pyplot as plt
import numpy as np
from math import exp # exponential function
dy = lambda x,y: x*y
f = lambda x: exp(x**2/2) # analytical solution function
x = 0
xn = 2
y = 1
# Container for step sizes
h = [0.5,0.2,0.1]
# Container to store the x values at each stepsize
# X =np.zeros((3,))
print ('x \t\t y (Euler) \t y (analytical)')
print ('%f \t %f \t %f'% (x,y,f(x)))
for j in range(0,len(h),1):
n = int((xn-x)/h[j])
for i in range(n):
y += dy(x, y)*h[j]
x += h[j]
print ('%f \t %f \t %f'% (x,y,f(x)))
plt.plot(x,y)
plt.show()
x y (Euler) y (analytical)
0.000000 1.000000 1.000000
0.500000 1.000000 1.133148
1.000000 1.250000 1.648721
1.500000 1.875000 3.080217
2.000000 3.281250 7.389056
所以问题实际上是尝试为不同的步长创建欧拉方法,即“如何更改我们的函数以循环遍历列表并使用 matplotlib 绘制结果”?
【问题讨论】:
您需要为 x 和 y 值创建一个列表,并在每一步附加当前值。当您为不同的 h 开始新的迭代时,不要忘记重置列表。 【参考方案1】:你犯了一个小错误,如果你想绘制结果,你需要将结果存储在一个容器中。我稍微重写了你的代码。在讨论您的代码有什么问题之前,我首先给您完整的代码。也许你自己发现了错误。我还添加了解析解的计算和其他一些你喜欢的小改进。所以这里是代码:
import matplotlib.pyplot as plt
import numpy as np
from math import exp # exponential function
dy = lambda x, y: x * y
f = lambda x: exp(x ** 2 / 2) # analytical solution function
x_final = 2
# analytical solution
x_a = np.arange(0, x_final, 0.01)
y_a = np.zeros(len(x_a))
for i in range(len(x_a)):
y_a[i] = f(x_a[i])
plt.plot(x_a, y_a, label="analytical")
# Container for step sizes
h = [0.5, 0.2, 0.1]
for j in range(len(h)):
x = 0
y = 1
print("h = " + str(h[j]))
print("x \t\t y (Euler) \t y (analytical)")
print("%f \t %f \t %f" % (x, y, f(x)))
n = int((x_final - x) / h[j])
x_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y
for i in range(n):
y += dy(x, y) * h[j]
x += h[j]
print("%f \t %f \t %f" % (x, y, f(x)))
x_n[i + 1] = x
y_n[i + 1] = y
plt.plot(x_n, y_n, "x-", label="h=" + str(h[j]))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
这会在我的计算机上绘制以下内容:
请注意,我将您的变量 xn
重命名为 x_final
以避免名称与我介绍的变量混淆。如前所述,您需要将每个 x 和 y 值存储在容器中。我为此使用了 NumPy 数组,但您也可以使用列表。这个
n = int((x_final - x) / h[j])
x_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y
只创建 2 个零数组,它们的大小等于子步数 +1。然后我将第一个值设置为等于初始值。这必须在 h
的循环内,因为每个 h 的子步数 n
是不同的。
在您的i
-loop 结束时,我只是将当前的x
和y
值写入数组中的正确位置。
for i in range(n):
y += dy(x, y) * h[j]
x += h[j]
print("%f \t %f \t %f" % (x, y, f(x)))
x_n[i + 1] = x
y_n[i + 1] = y
您需要将数组传递给函数:
plt.plot(x_n, y_n, "x-", label="h=" + str(h[j]))
我还添加了一个将显示在图例中的标签并将线型更改为"x-"
。
您犯的一个错误导致您的 i
循环仅针对第一个 h
执行是您没有将 x
和 y
重置为其初始值。所以你的n
在第一次运行外循环之后总是0
。
您当然可以优化一些东西,比如使用类似
的东西for h in h_list:
...
这将比总是使用h[j]
而不仅仅是h
更具可读性,但我认为现在这已经足够了。 ;)
【讨论】:
【参考方案2】:是的,你做对了,但它们是错误的。而不是 yn+1 y Euler 你已经为不同的 n 值打印了 yn。根据你的代码正在编写实际上是 yn 的 y 欧拉。 所以我这边的一个建议是你可以把这个表做成这样的 xn , yn ,slope , yn+1 然后是分析表。 让你的代码完美无缺
【讨论】:
感谢@RavneetKaur,继续努力:)以上是关于不同步长的欧拉方法。如何更改算法的代码以考虑不同的步长值?的主要内容,如果未能解决你的问题,请参考以下文章
如何在使用 LSTM 神经网络时考虑不同时间间隔的罕见事件?