2d wave equation with finite differences blowing up

Question

I am (naively) trying to solve the 2d wave equation with finite differences. But the system blows up instantly.

For simplicity I set the constant $c=1$, then I am left with $$\Delta u =u_{tt}.$$

I set up the FDM with

N = 9
u = np.zeros((N**2,1))
u.reshape(N,N)[N//2,N//2] = 1

Now, at the time t=0 the function looks like this

u_0 = 
[[0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]]

To which I apply the five point diffence stencil Lh

import scipy.sparse as sp
import numpy as np
h = 1/(N+1)
T = np.diag(N*[-2]) + np.diag((N-1)*[1],k=-1) +np.diag((N-1)*[1],k=1)    
spt = sp.coo_matrix(T)
spi = sp.eye(N)
Lh = (sp.kron(spi,spt) + sp.kron(spt,spi) ) / h**2

Since Lh is my discretized Laplacian, for every timestep I add the acceleration that it gives me to some book-keeping vector v to change the velocity accordingly:

v = np.zeros_like(u)
print(u.reshape(N,N).astype(int))
for i in range(3):
    v += Lh.dot(u)
    u += v
    print(u.reshape(N,N).astype(int))

Here is what happens:

[[   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0   99    0    0    0    0]
 [   0    0    0   99 -398   99    0    0    0]
 [   0    0    0    0   99    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]]

In the next timestep:

[[     0      0      0      0      0      0      0      0      0]
 [     0      0      0      0      0      0      0      0      0]
 [     0      0      0      0   9999      0      0      0      0]
 [     0      0      0  19999 -79699  19999      0      0      0]
 [     0      0   9999 -79699 198800 -79699   9999      0      0]
 [     0      0      0  19999 -79699  19999      0      0      0]
 [     0      0      0      0   9999      0      0      0      0]
 [     0      0      0      0      0      0      0      0      0]
 [     0      0      0      0      0      0      0      0      0]]

What is the problem here? I have not set up a boundary condition to my solution. Maybe I should?

greetings

Have you tried with smaller time step size dt? Then, the code reads Lh.dot(u)dt and vdt. — Hui Zhang, Dec 02 '18 at 14:07
The wave speed is 1. You should use 1, but better at least 2 time steps for any space step, thus $dt < dx/2$. — Lutz Lehmann, Dec 02 '18 at 16:15
To be at least qualitatively close the wave has to be able to travel at speed $c$. With the Euler method that you employ, changes propagate at speed $dx/dt$. Thus you need $dx>c,dt$, and the larger the gap the better the quantitative accuracy. With implicit methods the picture is not as clear, as the inverse of a sparse matrix is usually not sparse, and you would have to look at the eigenvalues of the propagation matrix. — Lutz Lehmann, Dec 03 '18 at 13:07

2d wave equation with finite differences blowing up

0 Answers0