Блог

так что я увидел этот код для RK4 в стеке и нашел его очень полезным. Однако я не могу придумать способ построения графика для каждого значения y при каждом приращении(h) x.

 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))
 

Вы можете добавить каждую точку в список в виде кортежа, а затем выполнить операцию построения линии в списке кортежей. Вы можете найти его в прокомментированном коде ниже.

 import matplotlib.pyplot as plt


def f(x, y):
    return 2 * x ** 2 - 4 * x   y


def RK4(x0, y0):
    pts = []  # empty list
    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
        pts.append((x0, y0))  # appending the tuple
    plt.plot(*zip(*pts))  # plotting the list of tuple
    plt.show()
    return y0


b = 3
h = 0.001
print(RK4(1, 0.7182818))
 

Вы можете увидеть сюжет следующим образом

С точки зрения проектирования было бы предпочтительнее, если бы код RK4 и код построения были разделены, численное решение не должно быть связано с тем, как его результаты будут использоваться впоследствии.

Тогда следующее решение будет касаться построения массива времени, он может быть передан методу RK4 или построен внутри и возвращен, оба имеют преимущества. Если вас беспокоит скорость, массивы должны быть построены явно в их окончательном виде (см. Пример на math.SE), для целесообразности можно также строить их постепенно. Таким образом, код может быть изменен следующим образом

 def RK4(f,x0,y0,xb,dx):
    x, y = [x0],[y0]
    while x0 < xb:
        k1 = dx*f(x0,y0)
        k2 = dx*f(x0 0.5*dx,y0 0.5*k1)
        k3 = dx*f(x0 0.5*dx,y0 0.5*k2)
        k4 = dx*f(x0 dx,y0 k3)
        y0  = (k1 2*k2 2*k3 k4)/6
        x0  = dx
        x.append(x0); y.append(y0) # for vector y use y0.copy()
    return x,y
 

а потом позвоните как

 x,y = RK4(f=f,x0=1.0,y0=0.7182818,xb=3.0,dx=1e-3)
plt.plot(x,y)
#title, axis labels
plt.grid(); plt.show()