import math as m
#import numpy as np

def myfunc(x):
    return m.cos(x)*x


def f(x, g):
    return m.cos(x) - m.sin(x)*x + k * (g - myfunc(x))


def absErr(arr):
    abserr = [abs(myfunc(x[i]) - arr[i]) for i in range(0, N+1)]
    return abserr


def globerr(arr):
    return max(absErr(arr))


def k1(i, t):
    return f(x[i-1], t[i-1])


def k2(i, t):
    return f(x[i-1] + h/2, t[i-1] + k1(i, t) * h/2)


def k3a(i, t):
    return f(x[i-1] + h, t[i-1] - h * k1(i, t) + 2 * h * k2(i, t))


def k3b(i, t):
    return f(x[i-1] + h / 2, t[i-1] + h * k2(i, t) / 2)


def k4(i, t):
    return f(x[i-1] + h, t[i-1] + h * k3b(i, t))


def euler():
    euler = start.copy()
    for i in range(1, N+1):
        euler[i] = euler[i-1] + h * f(x[i-1], euler[i-1])
    return euler


def rka():
    rka = start.copy()
    for i in range(1, N+1):
        rka[i] = rka[i-1] + h * (k1(i, rka) + 4 * k2(i, rka) * k3a(i, rka)) / 6
    return rka


def rkb():
    rkb = start.copy()
    for i in range(1, N+1):
        rkb[i] = rkb[i-1] + h * (k1(i, rkb) + 2 * k2(i, rkb) + 2 * k3b(i, rkb) + k4(i, rkb)) / 6
    return rkb


def adams_a():
    adams_a = startrka.copy()
    for i in range(4, N+1):
        adams_a[i] = adams_a[i-1] + (h/24) * (55 * f(x[i-1], adams_a[i-1]) - 59 * f(x[i-2], adams_a[i-2]) +
                                              37 * f(x[i-3], adams_a[i-3]) - 9 * f(x[i-4], adams_a[i-4]))
    return adams_a


def adams_b():
    adams_b = startrkb.copy()
    for i in range(4, N+1):
        adams_b[i] = adams_b[i-1] + (h/24) * (55 * f(x[i-1], adams_b[i-1]) - 59 * f(x[i-2], adams_b[i-2]) +
                                              37 * f(x[i-3], adams_b[i-3]) - 9 * f(x[i-4], adams_b[i-4]))
    return adams_b


def abm_a():
    predicta = adams_a().copy()
    abm_a = startrka.copy()
    for i in range(4, N+1):
        abm_a[i] = abm_a[i-1] + h * (9 * f(x[i], predicta[i]) + 19 * f(x[i-1], abm_a[i-1]) -
                                     5 * f(x[i-2], abm_a[i-2]) + f(x[i-3], abm_a[i-3]))
    return abm_a


def abm_b():
    predictb = adams_b().copy()
    abm_b = startrkb.copy()
    for i in range(4, N+1):
        abm_b[i] = abm_b[i-1] + h * (9 * f(x[i], predictb[i]) + 19 * f(x[i-1], abm_b[i-1]) -
                                     5 * f(x[i-2], abm_b[i-2]) + f(x[i-3], abm_b[i-3]))
    return abm_b


def gir_a():
    gir_a = startrka.copy()
    for i in range(4, N+1):
        gir_a[i] = (48 * gir_a[i-1] - 36 * gir_a[i-2] + 16 * gir_a[i-3] -
                    3 * gir_a[i-4] + 12 * h * m.cos(x[i]) - m.sin(x[i])*x[i] - k * myfunc(x[i])) / (25 - k)
    return gir_a


def gir_b():
    gir_b = startrkb.copy()
    for i in range(4, N+1):
        gir_b[i] = (48 * gir_b[i-1] - 36 * gir_b[i-2] + 16 * gir_b[i-3] -
                    3 * gir_b[i-4] + 12 * h * (m.cos(x[i]) - m.sin(x[i])*x[i] - k * myfunc(x[i]))) / (25 - k * 12 * h)
    return gir_b


def pred_corr():
    pred = start.copy()
    corr = [0 for i in range(0, N+1)]
    corr[0] = pred[0]
    stop = 1e-2
    for i in range(1, N+1):
        while abs(corr[i] - pred[i])/abs(pred[i]) <= stop:
            pred[i] = pred[i-1] + h * f(x[i-1], pred[i-1])
            corr[i] = corr[i-1] + (h/2) * (f(x[i-1], corr[i-1]) + f(x[i], pred[i]))


k = 5
N = 1000
h = 2 / N
x = [-1 + i * h for i in range(N+1)]
y = [0 for i in range(N+1)]
y[0] = myfunc(x[0])
start = [0 for i in range(0, N+1)]
start[0] = myfunc(x[0])
startrka = [0 for i in range(0, N+1)]
startrkb = [0 for i in range(0, N+1)]
for i in range(0, 4):
    startrka[i] = rka()[i]
    startrkb[i] = rkb()[i]
#plot = open('J:\Test2\ ' + str(N) + '.txt', 'w')
#for i in range(N+1):
#    plot.write(str(x[i]) + ' ' + str(myfunc(x[i])) + ' ' + str(adams_a()[i]) + '\n')
#plot.close()
test = open('J:\Test\ ' + str(N) + '.txt', 'w')
test.write(str(globerr(euler())) + ' ' + str(globerr(rka())) + ' ' + str(globerr(rkb())) + ' ' +
           str(globerr(adams_a())) + ' ' + str(globerr(adams_b())) + ' ' + str(globerr(gir_b())) + '\n')
test.close()