Periodic boundary condition and Runge-Kutta method of order 4 used in solving Partial Differential Equation(1D Wave packet equation).

415 Views Asked by At

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.