Euler explicit in Fourier space

Question

I am trying to implement an explicit Euler scheme in the Fourier space in order to solve the unsteady heat equation in one dimension :

The scheme is written as :

With the Fourier series of the function f written as :

I wrote k , the wave number, as a value that "moves" along the values of this list :

liste = list(range(int(N/2)+1))+list(range(-int(N/2)+1,0))

And it works very well (I have already tested it in previous simulations).

In order to implement the Euler scheme I have created this function that should build the values of the coefficient step time by step time :

def update(Time,dt) : 
    
    liste = list(range(int(N/2)+1))+list(range(-int(N/2)+1,0))
    nt = int(Time/dt)
    T = np.zeros(N)
    
    for i in range(N):
        T[i]  = coef_fourier_initial[i]
    
    for k in range(nt):
        for i in range(N):
            
            H = T[i]*(-4*(np.pi**2)*liste[i]**2)
            tmp   = T + dt*H 
            T     = tmp 
    
    return T

Unfortunately when I try to run it for :

y_2 = update(0.1,0.001)

I have these none values that appear and I don't understand why ?

[nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan

nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan

nan nan nan nan]

In order to simplify the comprehension of my code here I put how I've implemented the values of the function at t = 0 :

N = 40 
x = np.linspace(0,1-(1/N),N) 
l = np.arange(N)
f = np.zeros(N)

for i in l:
    if i>=0 and i<=9 :
        f[i] = 0
    elif i>=10 and i<=19 :
        f[i] = 2*(x[i] - 0.25)
    elif i>=20 and i<=29 :
        f[i] = 2*(-x[i] + 0.75)
    elif i>=30 and i<=39:
        f[i] = 0

Then the base matrix MSPPH :

def MSPPH(N):

    matrice=np.zeros((N,N),dtype="complex")
    k=0
    liste=list(range(int(N/2)+1))+list(range(-int(N/2)+1,0))
    liste_imp=list(range(int((N-1)/2)+1))+list(range(-int((N-1)/2),0))

    if N%2==0:
        for t in liste:
            for i in range(N):
                x=i/N
                matrice[k,i]=np.exp(-2*t*np.pi*(1j)*x)
            k+=1

    else:
        for t in liste_imp:
            for i in range(N-1):
                x=i/N
                matrice[k,i]=np.exp(-2*t*np.pi*(1j)*x)
            k+=1

    return np.array(matrice,dtype="complex")

And how I have built the initial coefficient :

coef_fourier_initial = np.zeros(N) 
coef_fourier_initial = np.dot(MSPPH(N),f)

Could someone help me to understand my error ?

Thanks a lot.

Euler explicit in Fourier space

0 Answers