Python Forum

Hello everyone,
I am in a MSc 1st year in physics and I have to realize a code in Python 3.8 on Spyder in order to solve the equation of Klein-Gordon. In order to simplifie the problem I consider the 2D case so I just take into account 1 dimension of space et 1 dimension of time. Here is the equation of Klein-Gordon I have to solve :

-ℏ²ðttu(x,t)=-c²ðxxu(x,t)-(m²c4/ℏ²)u(x,t) where ðtt is the second derivative in time
ðxx is the second derivative in space
u(x,t) is the wave equation
In order to code it in Python I am using the finite difference method and therefore I discretize the derivatives, with n the index for time and j the index for space, as follows :
ðttu(xj,tn)=ðttujn=(ujn+1-2ujn+ujn-1)/dt²
ðxxu(xj,tn)=ðxxujn=(uj-1n-2ujn+uj-1n)/dx²

The scheme obtained is the implicit Euler scheme which is expressed by :
ujn+1=-c²(dt²/dx²)(uj+1n-2ujn+uj-1n)-(m²c4dt²/ℏ²)ujn-ujn-1+2ujn

Therefore my code is the following :

import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d
from cmath import *
import matplotlib.colors as colors
from mpl_toolkits.mplot3d import Axes3D
from math import *

#Constant used
c = 3.0*10**8 #speed of light in the void in m.s-1
c2=c**2
c4=c2**2
c6=c2**3
hbar = 1.1*10**-34 #reduced Planck constant in J.s
hbar2=hbar**2
i = complex(0,1) #definition of the complex
m = 9.1*10**-31 #electron mass in kg
m2=m**2
m4=m2**2

#Mesh for time and space
Lx = 3.0 #total length
tf = 50 #final time
Nt = 250 #number of iteration in time
Nx = 250 #number of iteration in space

dt = tf/(Nt+1) #time step
dx = Lx/(Nx+1) #space step
dt2=dt**2
dx2=dx**2

#Figure
x_array=np.linspace(0.0,Lx,Nx)
t_array=np.linspace(0.0,tf,Nt)
X,T=np.meshgrid(x_array,t_array)


#Code
    #Initial conditions 
psi0=np.zeros((Nt,Nx))#Dirichlet conditions on the wave equation psi

    #Euler implicit 
for t in range (0,Nt-1):
    PSI=psi0
    PSI[0]=PSI[-1]=psi0[Nx-1] #Boundaries conditions
    for x in range (0,Nx-1):
        PSI[t][x+1]=PSI[t][x-1] #Boundaries conditions
        PSI[t+1][x]=-c2*(dt2/dx2)*(PSI[t][x+1]-2*PSI[t][x]+PSI[t][x-1])-(m2*c4*dt2/hbar2)*PSI[t][x]-PSI[t-1][x]+2*PSI[t][x]
    psi=PSI

#Spliting the real part and imaginary part of psi 
psire=psi.real
psiim=psi.imag

#Ploting the real part of psi
fig=plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('X', fontsize = 16)
ax.set_ylabel('time', fontsize = 16)
ax.set_zlabel('Wave function', fontsize = 16)
ax.view_init(elev=15, azim = 120)

p = ax.plot_surface(X,T,psire,cstride=1,linewidth=0,cmap='jet')

This code gives nothing but a blue square in 0 which is strange. I thought this might come from my initial and boundaries conditions or from the method I use or from the code itself but I tried several modifications (you can find them below) which just gives either nothing or another blue square. I am quite lost and I wonder if you could help me to solve this problem. Here are the modifications I made :
1) Initial condition as a wave packet without dispersion : the code for the wave packet only comes from https://physique.ensc-rennes.fr/tp_simulation_ondes.php

import numpy as np
from matplotlib import *
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d
from cmath import *
import matplotlib.colors as colors
from mpl_toolkits.mplot3d import Axes3D
from math import *

#Constant
c = 3.0*10**8 #speed of light in m.s-1
c2=c**2
c4=c2**2
hbar = 1.1*10**-34 #reduced Planck constant in J.s
hbar2=hbar**2
i = complex(0,1) #definition of the complex
m = 9.109*10**-31 #electron mass in kg
m2=m**2


#Definition od the mesh of space and time
Lx = 3.0 #total length
tf = 50 #total time
Nt = 250 #Number of time step
Nx = 250 #number of space step 

dt = tf/(Nt+1) #time step
dx = Lx/(Nx+1) #space step
dt2=dt**2
dx2=dx**2

#Defintiion of the position of calculus
positionsx=[j*Lx/Nx for j in range(Nx)]

temps=[j*tf/Nt for j in range(Nt)]

#period
T=0.33

#courbes
mescourbes=[[exp(-(t-x/c)**2/1.1)*sin(2*pi/T*(t-x/c)) for x in positionsx] for t in temps]

#Figure
x_array=np.linspace(0.0,Lx,Nx)
t_array=np.linspace(0.0,tf,Nt)
X,T=np.meshgrid(x_array,t_array)

#Code 
psi0=mescourbes*np.ones((Nt,Nx)) #Initial conditions

    #Euler implicit
PSI=psi0
for t in range (0,Nt-1):
    PSI[0]=PSI[-1]=psi0[Nx-1]#Boundaries conditions 
    for x in range (0,Nx-1):
        PSI[t+1][x]=-c2*(dt2/dx2)*(PSI[t][x+1]-2*PSI[t][x]+PSI[t][x-1])-(m2*c4*dt2/hbar2)*PSI[t][x]-PSI[t-1][x]+2*PSI[t][x]
    psi=PSI

psire=psi.real
psiim=psi.imag

#Ploting the real part of psi
fig=plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('X', fontsize = 16)
ax.set_ylabel('time', fontsize = 16)
ax.set_zlabel('FWave function', fontsize = 16)
ax.view_init(elev=15, azim = 120)

p = ax.plot_surface(X,T,psire,cstride=1,linewidth=0,cmap='jet')

I don't know what to do anymore to solve this equation, can you help me please ? Maybe should I change my boundaries conditions into Cauchy conditions but I do not know how to code it ? Do you have any website that could help me please ?

Thanks for your time,
best regards

tdub