# 2d wave equation with finite differences blowing up

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 v*dt. – Hui Zhang Dec 2 '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$. – LutzL Dec 2 '18 at 16:15
• thank you,Hui Zhang. @LutzL why is that? – dba Dec 3 '18 at 12:23
• 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. – LutzL Dec 3 '18 at 13:07