Thread: Numerical solution to hyperbolic PDE - grid leapfrog.

1. Numerical solution to hyperbolic PDE - grid leapfrog.

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!

2. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

Originally Posted by kornellster
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!
1. It appears you have assumed that $c_{s},$ the "isothermal velocity of sound", is equal to $-1/2.$ Is that correct?

2. If I look at Equation 2.66, and translate it a bit in terms of your code variables, I get the following:

$B\left(j,n+\frac{1}{2}\right)=B\left(j,n-\frac{1}{2}\right)-c_{s}\,\frac{\Delta t}{\Delta x}\left[A\left(j+\frac{1}{2},n\right)-A\left(j-\frac{1}{2},n\right)\right]$

$A\left(j+\frac{1}{2},n+1\right)=A\left(j+\frac{1}{ 2},n\right)-c_{s}\,\frac{\Delta t}{\Delta x}\left[B\left(j+1,n+\frac{1}{2}\right)-B\left(j,n+\frac{1}{2}\right)\right].$

Some pieces in these two equations that seem to be missing in your code are the following:

A. Where is $c_{s}?$
B. Where is $\Delta x?$
C. Have you thought through the fractions in the indexing carefully? I just wonder if there's some scaling that's off somewhere.

These are a few thoughts I had.

3. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

I have not wrote it explicitly, but in the picture with the graphs from the book, underneath it is stated that $\Delta x=c_{s}=1$. So I guess that doesn't matter.

Does code for the scheme look right?

PS. I'm not sure about the part where I
Code:
 A_array[0,j+1]:=1.0:
before, when I forgot about it, my result was totally wrong... Now I'm just naively setting it to 1 but it seems right since the waves don't get there by the time I finish the integration...

Yay!

I kind of got it. The main problem seems to be not knowing the initial pressure distribution function. Using a gaussian of the same height and width gives the same mid-grid oscillation as in the book, but overall it looks worse (becomes very steep and ugly). Do you have alternatives for sin and gaussian what would have a similar shape?

4. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

Originally Posted by kornellster

I have not wrote it explicitly, but in the picture with the graphs from the book, underneath it is stated that $\Delta x=c_{s}=1$. So I guess that doesn't matter.

Does code for the scheme look right?
I'm not sure where this line:

B_array[k,0]:= dt/2.0*(A_array[k+1,0]-A_array[k,0]):

came from. Can you justify that line? That's the line where I came to the conclusion that you thought that $c_{s}=-1/2.$

PS. I'm not sure about the part where I
Code:
 A_array[0,j+1]:=1.0:
before, when I forgot about it, my result was totally wrong... Now I'm just naively setting it to 1 but it seems right since the waves don't get there by the time I finish the integration...
What precisely are your boundary conditions? Those are typically given to you, especially in a numerical problem.

Yay!

I kind of got it. The main problem seems to be not knowing the initial pressure distribution function. Using a gaussian of the same height and width gives the same mid-grid oscillation as in the book, but overall it looks worse (becomes very steep and ugly). Do you have alternatives for sin and gaussian what would have a similar shape?
I'm of the opinion that you should have been given the initial/boundary conditions. So you shouldn't have to guess what they are.

5. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

No,boundary conditions were not given. All I know is that the is a perturbation in the pressure (the initial pressure distribution) and that at t=0 velocities are 0. I talked with my lecturer and he said i have to figure it out, that it is the whole point of these assignments

That line you wanted me to justify is the initial Euler backward step to find velocities at t=-1/2, such that at t=0, I would get zero velocities...

6. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

Originally Posted by kornellster
No,boundary conditions were not given. All I know is that the is a perturbation in the pressure (the initial pressure distribution) and that at t=0 velocities are 0. I talked with my lecturer and he said i have to figure it out, that it is the whole point of these assignments

That line you wanted me to justify is the initial Euler backward step to find velocities at t=-1/2, such that at t=0, I would get zero velocities...
So, what does the alternative

A_array[0,j+1]:=0.0:

give you?

I would double-check the line there that I wanted you to justify. Something looks a little strange (to me) about it.

7. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

It gives a 0 pressure at the boundary point of the grid which acts as a source of another wave and kills the whole calculation... The result then gets really strange...

8. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

OK, I think I understand this slightly better. My oscillations are different than in the screenshot because my initial perturbation was symmetric and a majority must have cancelled out. Making the function slightly asymmetric gives exactly the same results as suggested in the book.

Do you know where can I find a Runge Kutta scheme for this type of partial differential equation?

9. Re: Numerical solution to hyperbolic PDE - grid leapfrog.

Originally Posted by kornellster
OK, I think I understand this slightly better. My oscillations are different than in the screenshot because my initial perturbation was symmetric and a majority must have cancelled out. Making the function slightly asymmetric gives exactly the same results as suggested in the book.

Do you know where can I find a Runge Kutta scheme for this type of partial differential equation?