欧拉近似

Question

我正在尝试为微分方程的欧拉近似制作一个程序（我在一本书中遇到了一个问题告诉我）。我设法编写了适用于所有事物的东西（至少我已经测试过）但是由于某种原因，它在 step = 0,1 失败...... 代码如下：

from math import pow as p

y0: int = 3

h = input("Step size:")
h =float(h)
h1 = h
x_n=0
while x_n < (1):
    yn = y0 + h*(3*(p(x_n ,2))*(2-y0))
    x_n=h1
    h1=h1+h
    y0=yn
else:
    print(yn)

# h(1),f(1)=3
# h(0.1),f(1)=2,39279
# h(0.01),f(1)=2.37011
# h(0.001),f(1)=2.36810

它适用于 h=1、h=0.5、h=0.2、h=0.05、h=0.01、h=0.001 等等......由于某种原因，它只在 h=0.1 时失败。

有问题的微分方程是 y'=3(x^2)(2-y)

感谢您的帮助！

Answer 1

您的代码没有问题。 （您也可以使用** 代替pow；例如pow(2,3) = 2**3 = 8

Answer 2

（如果变量名有一个统一的结构会更好看）

如果你也打印出最后一步的时间

for h in [1.0,0.5,0.2,0.1,0.05]:
    x_n, y_0 = 0.0, 3.0
    while x_n < (1):
        y_n = y_0 + h*(3*(x_n**2)*(2-y_0))
        x_n = x_n + h
        y_0 = y_n
    else:
        print(f"h={h}: f({x_n:12.10f}) = {y_n:12.10f}")

然后你在结果中找到

h= 1.00000: f(1.0000000000) = 3.0000000000
h= 0.50000: f(1.0000000000) = 2.6250000000
h= 0.20000: f(1.0000000000) = 2.4261034230
h= 0.10000: f(1.1000000000) = 2.2749557398
h= 0.05000: f(1.0000000000) = 2.3795619396

确实在h=0.1 迭代只在x=1.1 停止。这是由于变量的浮点格式，靠近x=1 的点x_n 可能很快落在下面，因此执行另一个步骤。防止这种情况的一种方法是调整最后一步（另一种是使用“近距离”比较）

for h in [1.0,0.5,0.2,0.1,0.05]:
    x_n, y_0, x_f = 0.0, 3.0, 1.0
    y_n = y_0
    while x_n < x_f:
        if x_n + 1.001*h > x_f: h = x_f - x_n
        y_n = y_n + h*(3*(x_n**2)*(2-y_n))
        x_n = x_n + h
    else:
        print(f"h={h:8.5f}: f({x_n:12.10f}) = {y_n:12.10f}")

结果更一致

h= 1.00000: f(1.0000000000) = 3.0000000000
h= 0.50000: f(1.0000000000) = 2.6250000000
h= 0.20000: f(1.0000000000) = 2.4261034230
h= 0.10000: f(1.0000000000) = 2.3927939140
h= 0.05000: f(1.0000000000) = 2.3795619396

---

为了探索数字的浮点表示的影响，打印出浮点数的真正含义，

print(", ".join("%20.17f"%x for x in np.cumsum(10*[0.1])))

给予

0.10000000000000001,  0.20000000000000001,  0.30000000000000004,  
0.40000000000000002,  0.50000000000000000,  0.59999999999999998,  
0.69999999999999996,  0.79999999999999993,  0.89999999999999991,  
0.99999999999999989

这表明浮点运算后最后尾数位的舍入可能会产生相当违反直觉的结果