Рунге-Кутта 4 в питоне - PullRequest
Новая
       15

Рунге-Кутта 4 в питоне

0 голосов
/ 25 апреля 2018

У меня есть вопрос, в коде, для 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))

1 Ответ

0 голосов
/ 25 апреля 2018

Пожалуйста, напечатайте также последний x0, вы увидите, что итерация никогда не останавливается на b. Представление с плавающей запятой h и увеличение x0 на h редко бывают точными, так что итерация сама по себе не "достигнет" значения b, она остановится на чуть меньше b+h. 

 h           returned   x       returned   y
------------------------------------------------
0.1     (3.0000000000000018, 2.0855312227119374)
0.01    (3.0099999999999794, 2.1671997121516187)
0.001   (3.0009999999997796, 2.093630295729702 )
0.0001  (3.000000000002,     2.08553671291036  )
1e-05   (3.0000000000131024, 2.0855367130274134)

Вы можете исправить это, предварительно вычислив число шагов или исправив последний шаг 

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,6): h=10**-k; print h, RK4(1,0.7182818,b,h)

, который дает желаемые результаты 

 h        returned (x, y)
-----------------------------------
0.1     (3.0, 2.0855312227119236)
0.01    (3.0, 2.0855367122312707)
0.001   (3.0, 2.085536712901802 )
0.0001  (3.0, 2.0855367128941893)
1e-05   (3.0, 2.0855367129214737)

Как было объявлено, вы можете предварительно вычислить число 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)
...