Реализация алгоритмов Метод прогноза и коррекции и Рунге - Кутта -Фельберга

import math
import numpy as np
 
def f(x,y):
    return x*x-2*y
 
def rkf45(x0,y0,e,h,a,b):
    n = int((b-a)/h)
    t = np.arange(a,b,h)
    yi = y0
    i = 0
    for i in range(n):
        k0 = h*f(t[i],yi)
        k1 = h*f(t[i]+(2/3)*h,yi+(2/3)*k0)
        k2 = h*f(t[i]+(1/3)*h,yi+(1/12)*k0+(1/4)*k1)
        k3 = h*f(t[i]+(3/4)*h,yi+(69/128)*k0-(143/128)*k1+(135/64)*k2)
        k4 = h*f(t[i]+h,yi-(17/22)*k0+(27/4)*k1-(27/5)*k2+(16/15)*k3)
        k5 = h*f(t[i]+(5/6)*h,yi+(65/432)*k0-(5/16)*k1+(13/16)*k2+(4/27)*k3+(5/144)*k4)
        yi = yi + (1/9)*k0+(9/20)*k2+(16/45)*k3+(1/12)*k4
        r = (1/150)*k0-(3/100)*k2+(16/75)*k3+(1/20)*k4-(6/25)*k5
        #print(f"\tIteration: {i}\tx: {t[i]}\tyi: {yi}\tError: {r}\tStep: {h}\n")
        if abs(r) > e:
            h = h/2
            yi = 1
        elif abs(r) < e/20:
            h = h*2
            yi = 1
        elif e/20<=abs(r)<=e:
            print(f"\tIteration: {i}\tx: {t[i]}\tyi: {yi}\tError: {r}\tStep: {h}\n")
    return yi
 
x0 = 0#float(input("x0: "))
y0 = 1#float(input("y0: "))
h = 0.1#float(input("h: "))
a = 0#float(input("a: "))
b = 1#float(input("b: "))
e = 0.001#float(input("e: "))
 
n = int((b-a)/h)
 
print(rkf45(x0,y0,e,h,a,b))

import math
import numpy as np
import matplotlib.pyplot as plt
 
def f(x,y):
    return x*x-2*y
 
def rk4(x0,y0,h,a,b):
    t = x0
    y = y0
    n = int((b-a)/h)
    ypoints = []
    tpoints = np.arange(a,b,h)
    for i in range(4):
        k1 = h*f(t,y)
        k2 = h*f(t+h*0.5,y+h*k1*0.5)
        k3 = h*f(t+h*0.5,y+h*k2*0.5)
        k4 = h*f(t+h,y+h*k3)
        y += (k1+2*k2+2*k3+k4)/6
        ypoints.append(y)
    
    return ypoints
 
def rkf45(x0,y0,h,a,b,e):
    y = rk4(x0,y0,h,a,b)
    x = np.arange(a,b,h)
    tpoints =[]
    for i in range(len(y)):
        tpoints.append(x[i])
    p = []
    pi = 0
    yi = 1
    y_ = []
 
    for i in reversed(range(len(y))):
        pi = y[i] + (h/24)*(55*f(x[i],y[i]) - 59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])-9*f(x[i-3],y[i-3]))
        yi = y[i] + (h/24)*(9*f(x[i+1],pi)+19*f(x[i],y[i])-5*f(x[i-1],y[i-1])+f(x[i-2],y[i-2]))
        p.append(pi)
        y_.append(yi)
        print(f"\tIteration: {i}\tPredictor: {pi}\tCorrector: {yi}\tError: abs(yi-pi) : {abs(yi-pi)}\t\n")
        if abs(yi-pi) < e:
            print("\t\t\t\t\t\t\t\tEND")
            print(f"\t\t\t\t\t\t\t\tERROR: {abs(yi-pi)}")
            print(f"\t\t\t\t\t\t\t\tyn = {f(x[i],yi)}")
            break
    plt.ylabel('rkf45')
    plt.plot(tpoints, y, '-o')
    plt.grid(True)
    plt.xticks(x)
    plt.yticks(y_)
    plt.show()  
    return y_
 
 
x0 = 0#float(input("x0: "))
y0 = 1#float(input("y0: "))
h = 0.1#float(input("h: "))
a = 0#float(input("a: "))
b = 1#float(input("b: "))
e = 0.001#float(input("e: "))
print("\t\t\t\t\t\t\t\tSOLUTION\t")
print("----------------------------------------------------------------------------------------------------------------------------------------------")
print(f"\tyi: {rkf45(x0,y0,h,a,b,e)}")
print("---------