import numpy as np
import matplotlib.pyplot as plt

def legendre(k,t):
    m=len(t)
    pen=np.ones(m)
    ult=t
    if k == 0:
        val=pen
    if k == 1:
        val=ult
    if k > 1:
        for j in range(2,k+1):
            val = (2*j-1)/j *t*ult - (j-1)/j *pen
            pen=ult
            ult=val
    return val

def L2_data(n,f):
    # zusammengesetzte Simpson-Formel mit n**2 Intervallen
    N = n**2
    h= 2/N
    fl=np.zeros(N)
    fr=np.zeros(N)
    fm=np.zeros(N)
    pl=np.zeros(N)
    pm=np.zeros(N)
    pr=np.zeros(N)
    for j in range(N):
        pl[j]=(-1+j*h)
        fl[j]=f(-1+j*h)
        pr[j]=(-1+(j+1)*h)
        fr[j]=f(-1+(j+1)*h)
        pm[j]=(-1+j*h+h/2)
        fm[j]=f(-1+j*h+h/2)

    a=np.zeros(n+1)
    b=np.zeros(n+1)
    for j in range(n+1):
        a[j]=h*sum(1/6*legendre(j,pl)**2+2/3*legendre(j,pm)**2+1/6*legendre(j,pr)**2)
        b[j]=h*sum(1/6*fl*legendre(j,pl)+2/3*fm*legendre(j,pm)+1/6*fr*legendre(j,pr))
    return a, b

def L2bestapprox(a,b,t):
    n=len(a)
    m=len(t)
    pen=np.ones(m)
    ult=t
    val=b[0]/a[0]*pen + b[1]/a[1]* ult
    for j in range(2,n):
            tmp = (2*j-1)/j *(t*ult) - (j-1)/j *pen
            pen=ult
            ult=tmp
            val = val + b[j]/a[j] * tmp
    return val

# numerischer Test
def f(x):
    return 1/(1+25*x**2)

n=19
a,b = L2_data(n,f)
res=500 #Aufloesung des Graphen
y=np.zeros(res)
for j in range(res):
    y[j]=-1+2*j/res

I_L2best=L2bestapprox(a,b,y)

true_f=np.zeros(len(y))
for j in range(len(y)):
    true_f[j]=f(y[j])

plt.plot(y, I_L2best, label="L2-best P_"+str(n))
plt.plot(y, true_f, label="f")
plt.legend(loc='lower center')
plt.show()