Hi!

I'm working on solving partial differential equations numerically (writing codes). And its literally killing my soul.

I'm to work out a problem that can be found in a book - I include the screenshots so that you exactly know what's going on. I don't understand why i'm getting different results than suggested by the author.

I'm solving eq. (2.63) (see picture below), trying to use scheme (2.66)

Below is the solution to a small amplitude sound wave obtained by the author of the book. Basically, the hump in the middle is the initial pressure, the two smaller >1 humps are final pressure distributions and the bigger humps are arbitrarily normalized velocity distributions. (so starting with an initial perturbation we get two smaller waves traveling in both directions).

for s=1 i get almost the same thing (without the oscillation in velocity on the right hand side). And for s=1.0001 i also get pretty much the same thing (as for s=1)... I start getting similar oscillations as those in the book after setting s=1.001

I get a similar result but the main differences are: I need a slightly higher s for my solution to start visibly breaking up, I don't have that massive oscillation in the central region of the grid. For the solution using s=1 I don't get that subtel oscillation in velocity on the right side of the grid...

Here's my attempt at the solution: (maple code). For the initial pressure distribution I used a cut-off sine function. To find the velocities I used a backwards euler step such that I would get v=0 at t=0.

Please note that A is the array with the phi values (pressure) and B is the array with the psi values (velocities). The spatial grid has 402 points, time grid has 142 points.

Code:

A_array:= Array(0..401,0..141):
B_array:=Array(0..401,0..141):
dt:=1.001;
for n from 0 to 160 do A_array[n,0]:=1; end do:
for n from 161 to 240 do A_array[n,0]:= 1+0.065*sin(Pi/80*(n-160.5)); end do:
for n from 241 to 401 do A_array[n,0]:= 1; end do:
for k from 0 to 400 do
B_array[k,0]:= dt/2.0*(A_array[k+1,0]-A_array[k,0]):
end do:
for j from 0 to 140 do
A_array[0,j+1]:=1.0:
for k from 0 to 400 do
B_array[k,j+1]:= B_array[k,j]-dt*(A_array[k+1,j]-A_array[k,j]):
end do:
for s from 0 to 400 do
A_array[s+1,j+1]:= A_array[s+1,j]-dt*(B_array[s+1,j+1]-B_array[s,j+1]):
end do:
end do:

Whoever tells me where I'm wrong or whats causing the different results, I'll award you with a Nobell prize!