import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as lal

N = 20
a = 0
b = np.pi/2
h = (b - a)/N
x = np.linspace(a,b,N + 1)
t = np.linspace(a,b,N + 1)


def K(x,t):
    K = np.sin(x + t)
    return K

def f(x):
    f = x
    return f

f_i = f(x)


K_ij = np.zeros((N + 1,N + 1))

for i in range(0,N + 1):
    for j in range(0,N + 1):
        if j == 0 or j == N:
            K_ij[i,j] = h/2*K(x[i],t[j])
        else:
            K_ij[i,j] = h*K(x[i],t[j])


#allg: A*x = b hier: A*y = f wobei f = x und y = x bzw t also A*y = t


Em = np.zeros((N + 1,N + 1))

for i in range(0,N + 1):
    for j in range(0,N + 1):
        while j == i:
            Em[i,j] = 1
            j += 1


A_ij = Em - K_ij
y = lal.solve(A_ij,f_i)
print(y)

plt.figure
plt.plot(x,y,label = 'N = 20')
plt.legend(loc = 'upper left')
plt.xlabel('x')
plt.ylabel('$y(x_{i})$')
plt.title('Numerische Lösung')
plt.grid()

def y_analytic(x):
    y_analytic = ((np.pi**2 - 4)*x -2*(4 - 2*np.pi + np.pi**2)*np.cos(x) - 8*(np.pi - 1)*np.sin(x))/(np.pi**2 - 4)
    return y_analytic

#print(y_analytic(N))


eps = 10**-4
N = 1
mdiff = 10
while eps < mdiff:
    h2 = (b - a)/N
    x2 = np.linspace(a,b,N + 1)
    t2 = np.linspace(a,b,N + 1)
    f_i2 = f(x2)
    K_ij2 = np.zeros((N + 1,N + 1))
    for i in range(0,N + 1):
        for j in range(0,N + 1):
            if j == 0 or j == N:
                K_ij2[i,j] = h2/2*K(x2[i],t2[j])
            else:
                K_ij2[i,j] = h2*K(x2[i],t2[j])
    Em2 = np.zeros((N + 1,N + 1))
    for i in range(0,N + 1):
        for j in range(0,N + 1):
            while j == i:
                Em2[i,j] = 1
                j += 1
    A_ij2 = Em2 - K_ij2
    y2 = lal.solve(A_ij2,f_i2)
    mdiff = max(abs(y2 - y_analytic(x2)))
#    print(mdiff)
#    print(N)
    N += 1
print('')
print('Anzahl der Intervalle um einen maximalen Fehler von',eps,'zu machen:',N)