不同步长的欧拉法。如何更改算法代码以说明步长的不同值？

Question

我有一个解决数值微分方程问题的算法，称为欧拉法。本质上，欧拉方法近似于微分方程的解。我的函数适用于单个步长（值 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

所以问题实际上是尝试为不同的步长创建欧拉方法，即 "How does one change our function to loop over a list and plot the results with matplotlib"?

Answer 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 循环结束时，我只是将当前的 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

不是用 x 和 y 调用 plt.plot，因为它们是标量，所以只绘制一个点，您需要将数组传递给函数：

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 更易读，但我认为现在这已经足够了。 ;)

Answer 2

是的，你做对了，但他们是小故障。而不是 yn+1 y Euler 你已经为不同的 n.As 值打印了 yn 你的代码正在写 y eulers 这实际上是 yn。所以我这边的一个建议是，你可以像 xn , yn ,slope , yn+1 这样的 table 然后是分析的。休息你的代码是完美的

不同步长的欧拉法。如何更改算法代码以说明步长的不同值？

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?

python

integration

for-loop

matplotlib

numerical-methods