我有一个问题，在代码中，对于 h=0.1 显示一个小错误，即 h=0.01 和 h=0.001。我不明白为什么？，但是当 h=0.0001 时，误差再次减小。
谢谢！
def f(x,y):
return 2*x**2-4*x+y
def RK4(x0,y0):
while x0 < b:
k1 = h*f(x0,y0)
k2 = h*f(x0+0.5*h,y0+0.5*k1)
k3 = h*f(x0+0.5*h,y0+0.5*k2)
k4 = h*f(x0+h,y0+k3)
y0+=(k1+2*k2+2*k3+k4)/6
x0+=h
return y0
b=3
h=0.001
print(RK4(1,0.7182818))
错误分析
如果您还打印最后一个x0，那么您会看到迭代永远不会在正好b停止。h的浮点表示将偏离机器 epsilon 的某个部分。如果它稍大，则迭代将执行正确数量的循环。如果它较小，则迭代将执行比必要多一个循环，并在略小于b+h处停止。
此外，该测试问题的线性DE具有易于计算的精确解
y' - y = f'(x) - f(x), f(x) = -2*x^2
=> (y(x)-f(x))*exp(-x) = (y0-f(x0))*exp(-x0)
使溶液的流动函数为
def phi(x, x0,y0): return (y0+2*x0**2)*np.exp(x-x0)-2*x**2
为原始代码提供结果
exact solution: 2.085536712902183
h returned x returned y to exact at b to exact at ret. x
------------------------------------------------------------------------------------
0.1 : 3.00000000000000178 2.08553122271193736 -5.49019e-06 -5.49019e-06
0.01 : 3.00999999999997936 2.16719971215161866 0.081663 -6.90252e-10
0.001 : 3.00099999999977962 2.09363029572970216 0.00809358 -3.81029e-13
0.0001 : 3.00000000000200018 2.08553671291035991 8.17701e-12 -7.99849e-12
1e-05 : 3.00000000001310241 2.08553671302741339 1.2523e-10 1.92886e-11
在过冲情况下，误差明显是h的一小部分，而其他线显示预期的4阶收敛与浮点误差累积竞争。
校正变体
您可以通过先验计算步骤数或更正最后一步来纠正此问题
def f(x,y):
return 2*x**2-4*x+y
def RK4(x0,y0,xf,h):
while x0 < xf:
if x0+h > xf: h=xf-x0
k1 = h*f(x0,y0)
k2 = h*f(x0+0.5*h,y0+0.5*k1)
k3 = h*f(x0+0.5*h,y0+0.5*k2)
k4 = h*f(x0+h,y0+k3)
y0+=(k1+2*k2+2*k3+k4)/6
x0+=h
return x0,y0
b=3
for k in range (1,5):
h0=10**-k;
for h in [2*h0, h0, 0.5*h0]:
xf,yf = RK4(x0, y0,b,h);
print(f'{h:6} : {xf:.17f} {yf:.17f} {yf-phi(b,x0,y0):12.6g} {(yf-phi(b,x0,y0))/h**4:12.6g}')
这给出了预期的结果
h returned x y error error/h^4
--------------------------------------------------------------------------------
0.2 : 3.00000000000000000 2.08546790477393795 -6.88081e-05 -0.0430051
0.1 : 3.00000000000000000 2.08553122271192359 -5.49019e-06 -0.0549019
0.05 : 3.00000000000000000 2.08553632856235494 -3.8434e-07 -0.0614944
0.02 : 3.00000000000000000 2.08553670239502909 -1.05072e-08 -0.0656697
0.01 : 3.00000000000000000 2.08553671223127068 -6.70912e-10 -0.0670912
0.005 : 3.00000000000000000 2.08553671285973685 -4.2446e-11 -0.0679137
0.002 : 3.00000000000000000 2.08553671290148124 -7.01661e-13 -0.0438538
0.001 : 3.00000000000000000 2.08553671290180187 -3.81029e-13 -0.381029
0.0005 : 3.00000000000000000 2.08553671290032216 -1.86073e-12 -29.7717
0.0002 : 3.00000000000000000 2.08553671290186093 -3.21965e-13 -201.228
0.0001 : 3.00000000000000000 2.08553671289418929 -7.99361e-12 -79936.1
5e-05 : 3.00000000000000000 2.08553671292377762 2.15947e-11 3.45516e+06
计算的误差系数表明，对于0.005和0.1之间的h,4阶方法误差占主导地位，对于较大的步长，高阶误差项太大，对于较低的h，必要步骤的数量增加了太多，以至于浮点误差的累积主导了方法误差。
如前所述，您可以预先计算步数N并确保N*h=xf-x0.为此，请替换
while x0 < xf:
if x0+h > xf: h=xf-x0
跟
Dx = float(xf-x0); N = int(0.5+Dx/h); h = Dx/N
for _ in range(N):
您仍然可以观察到x0中浮点误差的累积，
0.1 (3.0000000000000018, 2.0855312227119374)
0.01 (2.9999999999999796, 2.0855367122311055)
0.001 (2.9999999999997797, 2.085536712900021 )
0.0001 (3.000000000002, 2.08553671291036 )
1e-05 (3.0000000000131024, 2.0855367130274134)