import numpy as np
import copy
from scipy import optimize
import scipy.linalg as sla
from sympy import *

y0 = 1/10
t0 = 1
T = 2
N = 10
step = (T-t0)/N
s1 = 0.20230885
s2 = 0.34733149
nodeT = np.zeros(N + 1)
for ind in range(0, N + 1):
    nodeT[ind] = t0 + step * ind

def ExatSolve():
    u = Function('u')
    t = symbols('t')
    ans = dsolve(u(t).diff(t)-((3*u(t))/(2*t)+(3*t*(u(t)**(1/3)))/2), u(t))
    # ans = dsolve(u(t).diff(t)-(2*t*(u(t)**3))/(1-(t*u(t))**2), u(t))
    pprint(ans)
    Constant = solve([ans.rhs.subs(t, t0)-y0])
    pprint(Constant[0])
    C1 = symbols('C1')
    ans = ans.subs(Constant[0])
    pprint(ans)
    arr = np.zeros(N+1)
    for ind in range(N+1):
        arr[ind] = list(solve(ans.subs(t, nodeT[ind]))[0].items())[0][1]
        print("В точке t = {0}, значение u({0}) = {1}.".format(nodeT[ind], arr[ind]))
    return arr
exat_list = ExatSolve()


nodeT = np.zeros(N + 1)
for ind in range(0, N + 1):
    nodeT[ind] = t0 + step * ind

def Function(t, u):
    return (3*u)/(2*t)+(3*t*(u**(1/3)))/2
    # return (2*t*(u**3))/(1-(t*u)**2)
def EulerObv():
    solution = np.zeros(N+1)
    solution[0] = y0
    for ind in range(0, N):
        solution[ind+1] = solution[ind]+step*Function(nodeT[ind], solution[ind])
    return solution


def EulerUnObv():
    solution = np.zeros(N+1)
    solution[0] = y0
    for ind in range(0, N):
        solution[ind+1] = optimize.root(lambda ans: -ans +  solution[ind]+step*Function(nodeT[ind+1], ans), solution[ind]).x[0]
    return solution

def RassingOrderFull(size):
    B = np.zeros(size*2)
    for ind in range(size*2):
        B[ind] = 1 / (1 + ind)
    def equations(p):
        A = np.array(list(p)[:size])
        alpha = np.array(list(p)[size:])
        ans = np.zeros(size*2)
        for ind in range(size*2):
            ans[ind] = (A*(alpha**ind)).sum()
        ans = ans-B
        return tuple(ans)


    tup = optimize.root(equations, np.zeros(size*2)).x
    A = tup[:size]
    alpha = tup[size:]
    print(A)
    print(alpha)

def RaisingOrder():
    # def NormFunction(X, u):
    #     return Function(X+1)
    # q = 2
    # alpha = np.array([0, 1/2])
    # A = np.array([0, 1])
    solution = np.zeros(N + 1)
    solution[0] = y0
    for ind in range(0, N):
        solution[ind + 1] = solution[ind] + Function(nodeT[ind]+step/2, solution[ind]+step*Function(nodeT[ind], solution[ind])/2) * step
    return solution


def Bucher1(x, y):
    s = 2
    C = np.array([0, 2/3])
    A = np.array([[0, 0], [2/3, 0]])
    B = np.array([1/4, 3/4])
    # s = 1
    # C = np.array([0])
    # A = np.array([[0]])
    # B = np.array([1])

    K = np.zeros(s)
    for ind in range(s):
        to_y = (K*A).sum()
        K[ind] = Function(x+step*C[ind], y + step*to_y)
    return (K*B).sum()

def RungeKutBuch():
    solution = np.zeros(N + 1)
    solution[0] = y0
    for ind in range(0, N):
        solution[ind + 1] = solution[ind] + step * Bucher1(nodeT[ind], solution[ind])
    return solution

def Adams3():
    solution = np.zeros(N + 1)
    solution[0] = y0
    solution[1] = s1
    solution[2] = s2

    for ind in range(0, N-1):
        value = solution[ind+1] + step * ((2/3)*Function(nodeT[ind+1], solution[ind+1]) - (1/12)*Function(nodeT[ind], solution[ind]))
        solution[ind + 2] = optimize.root(lambda ans:
                                          ans-step*5/12*Function(nodeT[ind+2], ans)-value,
                                          solution[ind+1]).x[0]
    return solution

print("Решение искалось в точках")
print(nodeT)
print()

print("Метод Адамса третьего порядка")
res = Adams3()
print(res)
print("Невязка c Точным решение")
print(abs(res-exat_list))
print()
print()
print("Неявный метод Эйлера")
eul = EulerUnObv()
print(eul)
print("Невязка c методом А")
print(abs(eul-res))
print("Невязка c Точным решение")
print(abs(eul-exat_list))
print()
print("Явный метод Эйлера")
eul = EulerObv()
print(eul)
print("Невязка c методом А")
print(abs(eul-res))
print("Невязка c Точным решение")
print(abs(eul-exat_list))
print()
print("Метод Рунге-Кутта")
ruk = RungeKutBuch()
print(ruk)
print("Невязка c методом А")
print(abs(ruk-res))
print("Невязка c Точным решение")
print(abs(ruk-exat_list))
print()
print()
print("Последовательного повышения точности")
ro = RaisingOrder()
print("Невязка c методом А")
print(abs(ro-res))
print("Невязка c Точным решение")
print(abs(ro-exat_list))
print()