不同步长大小的欧拉法 。如何改变算法的代码，以考虑步长大小的不同值？

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

所以问题其实是想创建不同步长大小的欧拉法，即 "如何改变我们的函数在列表上循环，用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 因为子步数 n 对每个h来说都是不同的。

在你的最后 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

而不是调用 plt.plot 与 x 和 y由于它们是sclars，你需要将数组传递给函数，所以只是绘制一个单点。

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但我认为现在已经足够了。)