I have to code for the problem of solving 1D wave packet propagation with equation given as: $$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$$ The initial solution is given as: u(x, 0) = $$ e^{(\alpha(x−x0 )^2)}cos(k_0(x − x_0))$$ with ${\alpha} = -5 $ and ${x_0} = 5.0$.Tangent hyperbolic function(double sided) is used as grid. I need to solve the above equation using Runge Kutta of order 4 for time discretisation with periodic boundary conditions. I know that while getting matrices for solving derivatives is correct and the periodic boundary conditions are applied properly. But during application of RK4 method I am bit unsure if the periodic conditions are applied properly.
Spatial discretisation algo have derivative terms considering three points on RHS : j-1, j and j+1
N is the number of nodes in the grid.
Derivative term is written as ddisp in the code.
It is the part of code where RK4 is used.
# RK4
# disp = u0 # initial condition array
flag = True
#h = 0
timeee = 0.000
timemax = 5001 #ms
while(timeee<timemax):
predisp = np.zeros(N+4)
ddisp = np.zeros(N+2) # derivative of displacement
k1 = np.zeros(N)
K2 = np.zeros(N)
K3 = np.zeros(N)
K4 = np.zeros(N)
# stage 1
for l in range (0,N):
predisp[l+2] = u0[l]
predisp[N+1] = predisp[3]
predisp[N+2] = predisp[3]
predisp[N+3] = predisp[4]
predisp[3] = predisp[N]
predisp[0] = predisp[N-1]
for p in range (0,N):
ddisp[p+1] = uj_prime[p]
ddisp[N+1] = ddisp[3]
ddisp[N] = ddisp[3]
ddisp[0] = ddisp[N+1]
for i in range(1,N+1):
k1[i-1] = -c*ddisp[i]
k1[N-1] = k1[0]
for y in range (0,N-1):
u0[y] = predisp[y+2]+ 0.5*(dt/a_sec)*k1[y]
u0[N-1] = u0[0]
#Stage 2
u_initial_colmat = np.transpose(u0)
uj_prime = C @ u_initial_colmat
for k in range(0,N):
ddisp[k+1] = uj_prime[k]
ddisp[N+1] = ddisp[3]
ddisp[N] = ddisp[3]
ddisp[0] = ddisp[N+1]
for i in range(1,N+1):
K2[i-1] = -c*ddisp[i]
K2[N-1] = K2[0]
for y in range(0,N-1):
u0[y] = predisp[y+2] + 0.5*(dt/a_sec)*K2[y]
u0[N-1] = u0[0]
#Stage 3
u_initial_colmat = np.transpose(u0)
uj_prime = C @ u_initial_colmat
for k in range(0,N):
ddisp[k+1] = uj_prime[k]
ddisp[N+1] = ddisp[3]
ddisp[N] = ddisp[3]
ddisp[0] = ddisp[N+1]
for i in range(1,N+1):
K3[i-1] = -c*ddisp[i]
K3[N-1] = K3[0]
for y in range(0,N-1):
u0[y] = predisp[y+2] + 1.0*(dt/a_sec)*K3[y]
u0[N-1] = u0[0]
#Stage 4
u_initial_colmat = np.transpose(u0)
uj_prime = C @ u_initial_colmat
for k in range(0,N):
ddisp[k+1] = uj_prime[k]
ddisp[N+1] = ddisp[3]
ddisp[N] = ddisp[3]
ddisp[0] = ddisp[N+1]
for i in range(1,N+1):
K4[i-1] = -c*ddisp[i]
K4[N-1] = K4[0]
for y in range(0,N-1):
u0[y] = predisp[y+2] + (1/6)*(dt/a_sec)*(k1[y] + 2*K2[y] + 2*K3[y] + K4[y])
u0[N-1] = u0[0]
milli_time = timeee/a_sec
uu = ((np.exp(alpha*((x-c*milli_time-x0)**2))*(np.cos(k0*(x-c*milli_time-x0))))) # Analytical Solution
timeee = timeee + dt
if (timeee >= timemax):
exit()
Any help would be highly appreciated.