Различные аппроксимации в средневзвешенном для моей реализации rk4 все одинаковы

Question

Итак, я пытаюсь написать реализацию метода ступеньки Кутты 4-го порядка для численного приближения дифференциального уравнения, как часть моего курса по математике, но также и для изучения некоторого программирования, но проблема в том, что он использует эти коэффициенты внутри каждого шага приближения. Итак, мы начинаем со значения x и значения y, и оно хочет найти значение y для более позднего значения x, и я задаю ему размер шага, обычно 0,1, и он увеличивается на 0,1 в x и дает мне новое приближение для значение y, и на каждом из этих шагов он выполняет 4 приближения, которые я называю k1,k2,k3 и k4. Итак, он берет взвешенное среднее из этих 4 приближений и дает мне окончательное приближение, которое должно быть достаточно точным. Проблема в том, что для каждого шага я получаю k1=k2=k3, но k4 отличается. Это не кажется правильным. Я проверяю это на функции, кормя ее f(x,y)=2x-3y+1, используя размер шага h, и мои k следующие:

k1=f(x,y)
k2=f(x+.5h,y+.5hk1)
k3=f(x+.5h,y+.5hk2)
k4=f(x+h,y+hk3)

Так что, просто исследуя это алгебраически, чтобы получить k1 ​​= k2, я получаю 9y=1-6x, но исходные значения, которыми я его кормлю, равны (x0,y0)=(1,5), и, очевидно, 45!=1-6

Так что все это кажется неправильным. Я не получаю правильный ответ в соответствии с учебником, и эти k не должны быть одинаковыми. Так или иначе, вот мой код. Операторы печати между k просто работают как тест, чтобы увидеть, что k фактически изменяются, когда мы перебираем n.

import numpy

def rk4(x0,y0,xf,h,f):
    y=[]
    x=numpy.linspace(x0,xf,(xf-x0)/h+1)
    y.insert(0,y0)
    for n in range(len(x)-1):
        k1=f(x[n],y[n])
        print k1
        k2=f(x[n]+(1/2)*h,y[n]+(1/2)*h*k1)
        print k2
        k3=f(x[n]+(1/2)*h,y[n]+(1/2)*h*k2)
        print k3
        k4=f(x[n]+h,y[n]+h*k3)
        print k4
        y.insert(n+1,y[n]+(h/6)*(k1+2*k2+2*k3+k4))
        print y[n]+(h/6)*(k1+2*k2+2*k3+k4)

    print x
    print y

def twoxminusthreeyplusone(x,y):
    return 2*x-3*y+1

rk4(1,5,1.5,0.1,twoxminusthreeyplusone)

Итак, я получаю этот вывод

/usr/bin/python2.7 /home/t/PycharmProjects/untitled/chunk.py
-12.0
-12.0
-12.0
-8.2
3.86333333333
-8.39
-8.39
-8.39
-5.673
3.06961666667
-5.80885
-5.80885
-5.80885
-3.866195
2.52110925
-3.96332775
-3.96332775
-3.96332775
-2.574329425
2.14792644708
-2.64377934125
-2.64377934125
-2.64377934125
-1.65064553888
1.900100743
[ 1.   1.1  1.2  1.3  1.4  1.5]
[5, 3.8633333333333333, 3.0696166666666667, 2.5211092500000003, 2.1479264470833335, 1.9001007429979169]

Process finished with exit code 0

Так что я не понимаю. То же самое, что k, кажется главной проблемой. Я пробовал это несколькими различными способами, как прежде, вместо того, чтобы просто определять k как функции x[n],y[n]. Каждый k начинался как пустой список, такой как y, и заполнял его после каждого шага для циклического цикла по n с использованием.insert(n,x) или, возможно,.insert(n-1,x) в зависимости от того, как я определил каждый шаг, но это казалось излишне сложным, и я думаю, что на самом деле возникла та же проблема

Answer 1

У вас типичная проблема, как правило, неожиданная.

Заменить все ваши факторы (1/2) от 0.5 (или же (1/2)*h от h/2, но избегайте деления, где это возможно). Целочисленное деление приводит к 0 и, следовательно, к наблюдаемому поведению.

---

Но обратите внимание, что с Python 3 ваш код будет работать, в свою очередь, целочисленное поведение должно быть принудительным.