2
$\begingroup$

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

$\endgroup$
  • 2
    $\begingroup$ Have you tried with smaller time step size dt? Then, the code reads Lh.dot(u)*dt and v*dt. $\endgroup$ – Hui Zhang Dec 2 '18 at 14:07
  • $\begingroup$ The wave speed is 1. You should use 1, but better at least 2 time steps for any space step, thus $dt < dx/2$. $\endgroup$ – LutzL Dec 2 '18 at 16:15
  • $\begingroup$ thank you,Hui Zhang. @LutzL why is that? $\endgroup$ – dba Dec 3 '18 at 12:23
  • $\begingroup$ 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. $\endgroup$ – LutzL Dec 3 '18 at 13:07

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.