import numpy as np
import matplotlib.pyplot as plt
import math

N = 129
err = []
def error(analytical,numerical):
    return np.sum(abs((numerical[(N/2)-2]-analytical[(N/2)-2])))

def sigma(t,D):
    return math.sqrt(2*D*t)

def analytical(t,x,D):
    firstPart = 1/(sigma(t,D)*math.sqrt(2*math.pi))
    secondPart = math.exp(-((x-L/2)**2)/(2*sigma(t,D)**2))
    return firstPart*secondPart


testNs = np.linspace(20,200,20)
for i in testNs:
    N = int(i)


    L = 2*math.pi
    D = 0.1
    h = L / (N-1)

    dt = 0.97*((h**2)/(2*D))

    vonnewman = h**2/2*D
    totaltime = 2

    U = np.zeros(((int(totaltime/dt)),N))

    U[0][len(U[0])/2] = 1/h

    U[:,0] = 0
    U[:,-1] = 0

    for i in range(1,int(U.shape[0])):
        for j in range(1,int(U.shape[1])-1):
            U[i][j] = U[i-1][j]+((D*dt)/h**2)*(U[i-1][j-1]+U[i-1][j+1]-2*U[i-1][j])

    A = np.zeros(((int(totaltime/dt)),N))
    for i in range(1,int(U.shape[0])):
        iC = i * dt
        for j in range(int(U.shape[1])-1):
            jC = j * h
            A[i][j] = analytical(iC,jC,D)


    #plt.xlabel('Time')
    #plt.ylabel('Concentration')
    #plt.plot(A[-1])
    #plt.plot(U[-1])
    #plt.show()

    err.append(error(A[-1],U[-1]))

plt.plot(testNs,err)
plt.show()