[Tutor] Solve wave equation

[Steven D'Aprano]

David Craig wrote:
> Hi,
> I am trying to write some code that will solve the 2D wave equation by 
> the finite difference method, but it just returns an array full of zeros 
> and NaN's. Not sure where I am going wrong, the code is attached so if 
> could someone point me in the right direction I'd appreciate this.

I'm not enough of a mathematician to tell what the solution of the 2D wave 
equation by the finite difference method should be. Are you sure that an array 
full of zeroes and NANs is not the right answer? Perhaps this is telling you 
that there is no solution?

A couple of other random stumbles in the dark:

> dt = 1**-4

Do you realise that 1**-4 == 1?


> # Set up pressure field.
> p=empty([nx,nz,nsteps])

This sets p to an array with random values (whatever rubbish just happens to 
be in memory at the time). You are expected to initialize the values yourself.

Since you go on to try to initialize to all zeroes, you should consider using 
numpy's zeroes() function instead of empty().


> for t in range(0,nsteps-1):
>     for z in range(0,nz-1):
>         for x in range(0,nx-1):
>             p[x,z,t] = 0

I'm pretty sure you end the loops one iteration too early here. Your array is 
indexed by

x = 0...nx-1 inclusive
z = 0...nz-1 inclusive
t = 0...nsteps-1 inclusive

(that is, the half-open intervals [0, nx) etc.) but you initialize only the 
values x = 0...nx-2 inclusive, etc. This is because Python's range() is 
already half-open on the right:

py> range(0, 10)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

By subtracting one, you miss an item:

py> range(0, 10-1)
[0, 1, 2, 3, 4, 5, 6, 7, 8]



-- 
Steven