Results 1 to 9 of 9

Math Help - Numerical solution to hyperbolic PDE - grid leapfrog.

  1. #1
    Newbie
    Joined
    Nov 2011
    Posts
    10

    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!
    Last edited by mr fantastic; November 24th 2011 at 11:51 AM. Reason: Title.
    Follow Math Help Forum on Facebook and Google+

  2. #2
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2

    Re: Numerical solution to hyperbolic PDE - grid leapfrog.

    Quote Originally Posted by kornellster View Post
    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.
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Newbie
    Joined
    Nov 2011
    Posts
    10

    Re: Numerical solution to hyperbolic PDE - grid leapfrog.

    Thank you for your time.

    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?
    Last edited by kornellster; November 24th 2011 at 01:08 PM.
    Follow Math Help Forum on Facebook and Google+

  4. #4
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2

    Re: Numerical solution to hyperbolic PDE - grid leapfrog.

    Quote Originally Posted by kornellster View Post
    Thank you for your time.

    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.
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Newbie
    Joined
    Nov 2011
    Posts
    10

    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...
    Follow Math Help Forum on Facebook and Google+

  6. #6
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2

    Re: Numerical solution to hyperbolic PDE - grid leapfrog.

    Quote Originally Posted by kornellster View Post
    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.
    Last edited by Ackbeet; November 24th 2011 at 02:39 PM.
    Follow Math Help Forum on Facebook and Google+

  7. #7
    Newbie
    Joined
    Nov 2011
    Posts
    10

    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...
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Newbie
    Joined
    Nov 2011
    Posts
    10

    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?
    Follow Math Help Forum on Facebook and Google+

  9. #9
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2

    Re: Numerical solution to hyperbolic PDE - grid leapfrog.

    Quote Originally Posted by kornellster View Post
    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?
    Glad your problems seem to have ironed themselves out.

    I have no idea about using RK methods on pde's. RK methods were originally designed for ode's; I'm not sure you can even use them on pde's, except for ode aspects of them.
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Replies: 0
    Last Post: November 29th 2011, 03:59 PM
  2. Numerical solution in mathematica
    Posted in the Math Software Forum
    Replies: 0
    Last Post: May 12th 2010, 10:22 AM
  3. Numerical solution of a first order ODE
    Posted in the Differential Equations Forum
    Replies: 1
    Last Post: May 1st 2010, 08:52 AM
  4. Help with numerical solution of PDE
    Posted in the Differential Equations Forum
    Replies: 3
    Last Post: March 20th 2010, 01:37 PM
  5. numerical solution
    Posted in the Calculus Forum
    Replies: 3
    Last Post: May 7th 2008, 04:30 AM

Search Tags


/mathhelpforum @mathhelpforum