0
$\begingroup$

Im trying to solve the following 2D wave equation: $$u_{tt} = u_{xx} + u_{yy}, \hspace{3mm} u(x,y,0) = \cos(4 \pi x) \sin(4 \pi y), \hspace{3mm} u_t(x,y,0) = 0$$ with the periodic boundary condition in both spatial directions $$u(x+1,y,t) = u(x,y,t), \hspace{3mm} u(x,y+1,t) = u(x,y,t)$$ which has analytical solution $$\cos(4 \pi x) \sin(4 \pi y) \cos(\sqrt{32}\pi t)$$

I used finite differences to solve the problem, with second order central differences in time, $x$ and $y$ from which I obtain an explicit iteration scheme. Details of the implementation can be viewed here under the "implementation" section: https://hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book008.html

My issue is that Im not happy with the numerical accuracy of the method, and the running time of the program is also very high, which means I can't choose a large number of points in $x$- and $y$-directions to get higher accuracy as the code would take too long to process. Below is my code in Python and a convergence plot:


class wave_eq_2D:
    def __init__(self,init_m,init_n,init_t):
        self.M = init_m
        self.N = init_n
        self.T = init_t

        self.x_grid = np.linspace(0, 1, self.M[0] + 2)
        self.y_grid = np.linspace(0,1,self.N[0] + 2)
        self.t_grid = np.linspace(0, 1, self.T + 2)

        self.solution = []
        self.rel_error = []



    def wave_solver2(self,init_c):

        M = self.M
        N = self.N
        T = self.T

        h_x = self.x_grid[1] - self.x_grid[0]
        h_y = self.y_grid[1] - self.y_grid[0]
        k = self.t_grid[1] - self.t_grid[0]

        r_x = k / h_x
        r_y = k / h_y

        u = np.zeros(shape=(M + 2, N + 2, T + 2)) #initialize numerical solution u

        u[:, :, 0] = init_c(self.x_grid, self.y_grid) #apply initial condition g(x,y) = 4cos(4 * pi * x) * 4*sin(4* pi * x)

        for i in range(T+1): #iterate in time

            u_xx = u[2:, 1:-1, i] - 2 * u[1:-1, 1:-1, i] + u[:-2, 1:-1, i]
            u_yy = u[1:-1,2:,i] -2*u[1:-1,1:-1,i] + u[1:-1,:-2,i]
            if (i == 0):  # using Neumann condition at first time step

                u[1:-1,1:-1, i + 1] = u[1:-1,1:-1, i] + 0.5*(r_x ** 2) * u_xx + 0.5*(r_y**2)*u_yy #apply time scheme for n = 0

                #Apply periodic BCS using time scheme in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
                u[0, 1:-1, i + 1] = u[0,1:-1, i] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) + (r_y**2)/2*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
                u[-1,:,i+1] = u[0,:,i+1]
                u[1:-1, 0, i + 1] = u[1:-1, 0, i] + 0.5 * (r_x ** 2)/2 * (u[2:, 0, i] - 2 * u[1:-1, 0, i] + u[:-2, 0, i]) + (r_y ** 2)/2 * (u[1:-1, 1, i] - 2 * u[1:-1, 0, i] + u[1:-1, -2, i])
                u[:, -1, i + 1] = u[:, 0, i + 1]

            else:

                u[1:-1,1:-1, i + 1] = 2*u[1:-1,1:-1,i] -u[1:-1,1:-1,i-1] + (r_x ** 2)/2 * u_xx + (r_y**2)/2*u_yy #apply time scheme for general n

                # Apply periodic BCS using time scheme in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
                u[0, 1:-1, i + 1] = 2*u[0,1:-1,i] -u[0,1:-1,i-1] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) +(0.5*r_y**2)*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
                u[-1, :, i + 1] = u[0, :, i + 1]  # on the right side
                u[1:-1, 0, i + 1] = 2*u[1:-1,0,i] -u[1:-1,0,i-1] + (r_x ** 2) * (u[2:,0, i] - 2 * u[1:-1,0, i] + u[:-2,0, i]) +(r_y**2)*(u[1:-1,1,i] -2*u[1:-1,0,i] + u[1:-1,-2,i])
                u[:,-1,i+1] = u[:,0,i+1]


        self.solution = u[:,:,-1]

relative l2 error

The convergence plot also looks weird as the convergence rate is not smooth. This makes me think there is an error somewhere. Convergence rate is approximately of second order which is as expected.

Would appreciate any feedback on how I could improve my wave equation solver implementation.

Edit: This is my full code and how I test convergence by running several iterations and doubling the number of points in $x$- and $y$ each time:


import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

from plotting_tools import set_ax


class wave_eq_2D:
    def __init__(self,init_m,init_n,init_t):
        self.M = [init_m]
        self.N = [init_n]
        self.T = init_t

        self.x_grid = [np.linspace(0, 1, self.M[0] + 2)]
        self.y_grid = [np.linspace(0,1,self.N[0] + 2)]
        self.t_grid = np.linspace(0, 1, self.T + 2)

        self.solution = []
        self.rel_error = []



    def wave_solver2(self,init_c,iteration):

        M = self.M[iteration]
        N = self.N[iteration]
        T = self.T

        h_x = self.x_grid[iteration][1] - self.x_grid[iteration][0]
        h_y = self.y_grid[iteration][1] - self.y_grid[iteration][0]
        k = self.t_grid[1] - self.t_grid[0]

        r_x = k / h_x
        r_y = k / h_y

        u = np.zeros(shape=(M + 2, N + 2, T + 2))

        u[:, :, 0] = init_c(self.x_grid[iteration], self.y_grid[iteration])

        for i in range(T+1):

            u_xx = u[2:, 1:-1, i] - 2 * u[1:-1, 1:-1, i] + u[:-2, 1:-1, i]
            u_yy = u[1:-1,2:,i] -2*u[1:-1,1:-1,i] + u[1:-1,:-2,i]
            if (i == 0):  # using Neumann condition at first time step

                u[1:-1,1:-1, i + 1] = u[1:-1,1:-1, i] + 0.5*(r_x ** 2) * u_xx + 0.5*(r_y**2)*u_yy

                #Apply periodic BCS in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
                u[0, 1:-1, i + 1] = u[0,1:-1, i] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) + (r_y**2)/2*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
                u[-1,:,i+1] = u[0,:,i+1]
                u[1:-1, 0, i + 1] = u[1:-1, 0, i] + 0.5 * (r_x ** 2)/2 * (u[2:, 0, i] - 2 * u[1:-1, 0, i] + u[:-2, 0, i]) + (r_y ** 2)/2 * (u[1:-1, 1, i] - 2 * u[1:-1, 0, i] + u[1:-1, -2, i])
                u[:, -1, i + 1] = u[:, 0, i + 1]

            else:

                u[1:-1,1:-1, i + 1] = 2*u[1:-1,1:-1,i] -u[1:-1,1:-1,i-1] + (r_x ** 2)/2 * u_xx + (r_y**2)/2*u_yy

                # Apply periodic BCS in both spatial dimensions: u(x+1,y,t) = u(x,y,t) = u(x,y+1,t)
                u[0, 1:-1, i + 1] = 2*u[0,1:-1,i] -u[0,1:-1,i-1] + 0.5*(r_x ** 2) * (u[1,1:-1, i] - 2 * u[0,1:-1, i] + u[-2,1:-1, i]) +(0.5*r_y**2)*(u[0,2:,i] -2*u[0,1:-1,i] + u[0,:-2,i])
                u[-1, :, i + 1] = u[0, :, i + 1]  # on the right side
                u[1:-1, 0, i + 1] = 2*u[1:-1,0,i] -u[1:-1,0,i-1] + (r_x ** 2) * (u[2:,0, i] - 2 * u[1:-1,0, i] + u[:-2,0, i]) +(r_y**2)*(u[1:-1,1,i] -2*u[1:-1,0,i] + u[1:-1,-2,i])
                u[:,-1,i+1] = u[:,0,i+1]


        self.solution.append(u[:,:,-1])

        print("Solution", iteration + 1, "found.")




    def get_error(self, u, iteration):
        # Defining the discrete l2 norm.
        def disc_norm(v):
            return np.sqrt(np.sum(v ** 2) / len(v))

        # Defining a meshgrid from our x- and y-grids.
        x_disc = self.x_grid[iteration]
        y_disc = self.y_grid[iteration]
        X_disc, Y_disc = np.meshgrid(x_disc, y_disc)

        # Defining the numerical and analytical solution arrays
        U_disc = self.solution[iteration]
        u_disc = u(x_disc, y_disc,1)

        # Finding the relative error using the l2 norm.
        error_disc = disc_norm(u_disc - U_disc) / disc_norm(u_disc)
        # Saving the error data in the relative error list.
        self.rel_error.append(error_disc)

    def check_convergence(self, iterations, u, init_c):

        # Finding the l2 error of the initial system.
        self.wave_solver2(init_c, 0)
        self.get_error(u, 0)

        # Calculating the l2 errors of a few more iterations of the system.
        for i in range(1, iterations):
            # Doubling the size of the x-grid and y-grid.
            self.M.append(2 * self.M[i - 1])
            M = self.M[i]
            self.x_grid.append(np.linspace(0, 1, M + 2))
            self.N.append(2*self.N[i-1])
            N = self.N[i]
            self.y_grid.append(np.linspace(0,1,N + 2))

            # Solving for the new system.
            self.wave_solver2(init_c, i)
            self.get_error(u, i)

    def plot_convergence(self, line_order=2, space='x'):
        disc_error = self.rel_error
        print('M-values:', self.M)

        if space == 'x':
            gridsize = self.M
            x_label = r'$M$'
            line_label = r'$O(h^{%s})$'
        elif space == 'y':
            gridsize = self.N
            x_label = r'$N$'
            line_label = r'$O(k^{%s})$'
        else:
            return

        x_array = np.linspace(gridsize[0], gridsize[-1], 1000)

        def line(x, order, init_error):
            return init_error * x ** (-order) / (x[0]) ** (-order)

        print('M-values:', gridsize)
        print('disc error:', disc_error)

        fig = plt.figure()
        ax = fig.add_subplot(111)
        set_ax(ax,
               title='Convergence testing for 2D Wave Equation',
               xscale='log', yscale='log',
               x_label=x_label, y_label='Relative error')

        ax.plot(gridsize, disc_error, marker='o', label=r'$l_2$ norm error', lw=3, c='#1F77B4')
        ax.plot(x_array, line(x_array, line_order, disc_error[0]), label=line_label % line_order,
                lw=3, ls='--', c='#1F77B4')

        ax.legend(fontsize=15)
        plt.show()
        plt.close()
        

#Boundary condition g(x,y)
def g(x,y):
    return np.cos(4* np.pi * x) * np.sin(4 * np.pi * y)

#Analytical solution
def u(x,y,t):
    return np.cos(4*np.pi*x)*np.sin(4*np.pi*y)*np.cos(np.sqrt(32)*np.pi*t)

System = wave_eq_2D(70,70,1550)


System.check_convergence(4,u,g)
System.plot_convergence()



```
Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ I do not understand why you have computations in the boundary conditions. Shouldn't you just copy the already existing data to have a periodic border for the next Laplace computation, u[0,1:-1,i+1] = u[-2,1:-1,i+1]; u[-1,1:-1,i+1]=u[1,1:-1,i+1]; u[:,0,i+1]=u[:,-2,i+1]; u[:,-1,i+1]=u[:,1,i+1]? $\endgroup$
    – Lutz Lehmann
    Mar 22, 2021 at 20:38
  • $\begingroup$ You seem to mix 2 ideas for the periodic continuation. You use the identification of index 0 with index -1. For a grid with $M\times N$ rectangles, this would require that the subdivisions be constructed using length $M+1$ and $N+1$ in linspace. But you use $M+2$ and $N+2$, which would make sense if, as I assumed in my previous comment, you identify the first two indices [0,1] with the last two [-2,-1], so that $x_N=1$, $x_{N+1}=1+h_x$ etc. Then you also do not need the extra computations for row and column 0. $\endgroup$
    – Lutz Lehmann
    Mar 23, 2021 at 8:26
  • $\begingroup$ @LutzLehmann My idea was to have $x_0 = 0, x_{M+1} = 1$, but wouldn't then $x_0$ and $x_{M+1}$ be equal? With $x_0$ = $u[0]$ and $x_{M+1}$ being $u[-1]$ i.e point number $M+2$. $\endgroup$
    – Pame
    Mar 23, 2021 at 10:26
  • 1
    $\begingroup$ There is nothing wrong in the code in that regard, I just wonder about the semantics of $M$ and $N$. So if $M=N=10$, you get a grid of $11\times 11$ squares. Is that your intention, or is the $M+2$ a left-over from a Matlab implementation where the array indices start at $1$? $\endgroup$
    – Lutz Lehmann
    Mar 23, 2021 at 10:50

1 Answer 1

Reset to default
3
$\begingroup$

Regarding performance, Python is definitely the bottleneck. I have experienced the same issue with a 2D Euler code I had developed, even with vectorised operations everywhere possible. It was actually even worse, as I was using solve_ivp time schemes which reallocated memory at every step... You can try and profile your code to see where the bottlenecks are. A profiler is included in the Spyder IDE (shortcut F10). Can you share your code (with the convergence analysis) ?

Also, I see that there are sometimes (r_x ** 2)/2, sometimes only (r_x ** 2)in your code. Have you double-checked that this is correct? I remember that there was an issue with that in your previous question for the 1D case.

Regarding the convergence rate, maybe the grid spacings you tried are not fine enough, such that the true asymptotic convergence is not yet reached. Or maybe the error is not computed correctly. Also, the time step may affect the rate of convergence if it is not fine enough (such that the spatial error does not dominate).

I'll be happy to complement my answer if you give more information !

Improve this answer
$\endgroup$
1
  • $\begingroup$ I edited my post now to supply the full code where I calculate error and convergence plots. Regarding $r_x^2$ vs $r_x^2$/2 it is because the formula changes for $n = 0$. It is explained in the link I posted. Since the code runs slowly when choosing a high number of points in $x$ and $y$ means very high running time which is why I have to choose the low for the convergence plot. Im also wondering if you are able to achieve better results with your own revised version of this code? (without using temporal integration). $\endgroup$
    – Pame
    Mar 23, 2021 at 9:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.