Page 1 of 2 12 LastLast
Results 1 to 15 of 17

Math Help - 1D wave equation open boundary

  1. #1
    Newbie
    Joined
    Jul 2010
    Posts
    8

    1D wave equation open boundary

    I am trying to write a solver for a 1D wave equation in MATLAB, and I have run into interesting problem that I just can't find a way out of.

    I start with the wave equation, and then discretize it, to arrive at the following,
    <br />
U(n+1,j)=a*(U(n,j+1)-2*U(n,j)+U(n,j-1))+2*U(n,j)-U(n-1,j) (for (j=1...end-1))

    I'm trying to simulate an open end string (perturbate it in the middle for example, and I want the wave to dissapear on the border)

    Usually it is done (correct me if I'm wrong) with applying zero gradient boundary condition, which in my case is:
    U(n,end+1)=U(n,end-1)
    and leads to:
    U(n+1,end)=a*(-2*U(n,end)+2*U(n,end-1))+2*U(n,end)-U(n-1,end);

    Unfortunately this condition creates a partial reflection from the end which is supposed to be free.

    Please help

    Regards
    Last edited by emirs; July 8th 2010 at 12:14 AM.
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Member oldguynewstudent's Avatar
    Joined
    Oct 2009
    From
    St. Louis Area
    Posts
    244
    Quote Originally Posted by emirs View Post
    I am trying to write a solver for a 1D wave equation in MATLAB, and I have run into interesting problem that I just can't find a way out of.

    I start with the wave equation, and then discretize it, to arrive at the following,

    U{n+1}(j)=a*(U{n}(j+1)-2*U{n}(j)+U{n}(j-1))+2*U{n}(j)-U{n-1}(j) (for (j=1...end-1))

    I'm trying to simulate an open end string (perturbate it in the middle for example, and I want the wave to dissapear on the border)

    Usually it is done (correct me if I'm wrong) with applying zero gradient boundary condition, which in my case is:
    U{n}(end+1)=U{n}(end-1)
    and leads to:
    U{n+1}(end)=a*(-2*U{n}(end)+2*U{n}(end-1))+2*U{n}(end)-U{n-1}(end);

    Unfortunately this condition creates a partial reflection from the end which is supposed to be free.

    Please help

    Regards
    Can't read your code. Use TEX to put math tags around your items.

    U{n+1}(end)=a*(-2*U{n}(end)+2*U{n}(end-1))+2*U{n}(end)-U{n-1}(end);

    Then click preview post to see how it will look.
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Newbie
    Joined
    Jul 2010
    Posts
    8
    Is it OK now? Does anyone knows the answer?
    Follow Math Help Forum on Facebook and Google+

  4. #4
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    5
    Awards
    2
    I can't pretend to be an expert, but I'd that it's possible you're simply getting numerical error from the discretization. Are the reflections you're getting nearly 100% of the incident waves, or are they more like 1%?
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Member oldguynewstudent's Avatar
    Joined
    Oct 2009
    From
    St. Louis Area
    Posts
    244
    I'm not an expert on this equation, but I think Ackbeet is probably correct. I used to be an expert programmer before my job was shipped overseas so try these debugging tricks:

    1) check for rounding, to what accuracy do you have Matlab set?

    2) put displays in your code at the boundary conditions where you are getting the partial reflections like this: if (n=10, or 20, or j-10) then display certain variables.

    3) Matlab has excellent debugging features so try setting your breakpoints properly and step through the code line by line ( I know it's tedious, but sometimes necessary).

    4) try a simpler equation by copying your M-file and saving a work version you can mess with. In other words, the wave equation reduces to simpler versions for certain conditions so play with limiting factors to see if it behaves differently.

    5) Most important, go over your code and pretend you are explaining it to someone who doesn't know what you are doing. Do this line by line. I have found typos many times with this technique because your mind assumes you have written what you were thinking.

    6) If none of these work, attach your M-file to your post along with any data and I'll try to run it a couple of times to debug it.
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Newbie
    Joined
    Jul 2010
    Posts
    8
    I think the rounding is not an issue. Matlab works with double precision floating points, which I think is more than enough. The numerical method (using finite differences) is proven to be numerically stable, so it shouldn't give such results only because of rounding.
    Also the reflection of the waves is around 70% or more (cannot tell you exactly), but the reflection is not the usual solid boundary reflection (which causes the wave to change phase), but instead the waves look like remaining in phase.

    The code should be just fine. I will try to attach it.

    My guess that this has to have something to do with gradient and its numerical approximation, the speed of the waves (propagation speed) and the numerical discretization. The gradient on boundary shouldn't be zero, but something else. The propagation of the wave on x-axis is in some kind of relationship with the speed of oscillation around equilibrium. But I cannot find out what it is. I tried to use this code:
    U(n+1,end)=U(n,end-1);
    (what was before in end-1, now is in end), and it is giving a better result, but not a 100% better. It is missing something.
    Follow Math Help Forum on Facebook and Google+

  7. #7
    Newbie
    Joined
    Jul 2010
    Posts
    8
    wave1D.txt

    I have attached the code. Just rename .txt to .m file.
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Member oldguynewstudent's Avatar
    Joined
    Oct 2009
    From
    St. Louis Area
    Posts
    244
    Quote Originally Posted by emirs View Post
    wave1D.txt

    I have attached the code. Just rename .txt to .m file.
    I've run several simulations and stepped through the code.

    I have a question about line 38

    U{n+1}(end)=a*(-2*U{n}(end)+2*U{n}(end-1))+2*U{n}(end)-U{n-1}(end);

    When you have the term 2*U{n}(end-1) isn't that reflecting the wave back?

    If you want no reflection, wouldn't this just be U{n}(end-1)? Multiplying by two averages in the reflecting part doesn't it?

    Anyway I found the reflection to start a little before i=722.

    Try putting your breakpoint at line 45 drawnow.

    Hit the continue about 690-700 times until your wave gets near the boundary.

    Then double click on the U in the workspace, then double click on each U1, U2, and U3. Examine the values around columns 990-1001. Slowly step through each iteration until you see the end columns on U3 start to increase instead of decrease. Do the calculations in your formulas by hand at this point and the reason for the increase should become clear.

    Programming something complicated takes a lot of patience and intuition. Troubleshooting takes the highest skill and patience. I am sorry that I am bitter about loosing my job after 35 years, but I started out flipping binary switches and coding in octal. Now I'm thrown out like yesterday's garbage.

    Good luck!
    Follow Math Help Forum on Facebook and Google+

  9. #9
    Newbie
    Joined
    Jul 2010
    Posts
    8
    Quote Originally Posted by oldguynewstudent View Post
    I've run several simulations and stepped through the code.

    I have a question about line 38

    U{n+1}(end)=a*(-2*U{n}(end)+2*U{n}(end-1))+2*U{n}(end)-U{n-1}(end);

    When you have the term 2*U{n}(end-1) isn't that reflecting the wave back?

    If you want no reflection, wouldn't this just be U{n}(end-1)? Multiplying by two averages in the reflecting part doesn't it?

    Anyway I found the reflection to start a little before i=722.

    Try putting your breakpoint at line 45 drawnow.

    Hit the continue about 690-700 times until your wave gets near the boundary.

    Then double click on the U in the workspace, then double click on each U1, U2, and U3. Examine the values around columns 990-1001. Slowly step through each iteration until you see the end columns on U3 start to increase instead of decrease. Do the calculations in your formulas by hand at this point and the reason for the increase should become clear.

    Programming something complicated takes a lot of patience and intuition. Troubleshooting takes the highest skill and patience. I am sorry that I am bitter about loosing my job after 35 years, but I started out flipping binary switches and coding in octal. Now I'm thrown out like yesterday's garbage.

    Good luck!
    My english is unfortunately not so good so I couldn't understand what is your point.
    1.If I remove the two in 2*U{n}(end-1) part, then the code still doesn't work, and I still don't see any mathematical approval for that (I had shown where from that relation came in).
    2.I debbuged as you instructed me, but I still don't get your point?

    Thanks for the efforts
    Follow Math Help Forum on Facebook and Google+

  10. #10
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by emirs View Post
    I am trying to write a solver for a 1D wave equation in MATLAB, and I have run into interesting problem that I just can't find a way out of.

    I start with the wave equation, and then discretize it, to arrive at the following,
    <br />
U(n+1,j)=a*(U(n,j+1)-2*U(n,j)+U(n,j-1))+2*U(n,j)-U(n-1,j) (for (j=1...end-1))

    I'm trying to simulate an open end string (perturbate it in the middle for example, and I want the wave to dissapear on the border)

    Usually it is done (correct me if I'm wrong) with applying zero gradient boundary condition, which in my case is:
    U(n,end+1)=U(n,end-1)
    and leads to:
    U(n+1,end)=a*(-2*U(n,end)+2*U(n,end-1))+2*U(n,end)-U(n-1,end);

    Unfortunately this condition creates a partial reflection from the end which is supposed to be free.

    Please help

    Regards
    I think your condition on the supposedly free end is the problem, this looks like the you are constraining the rate of change of U at the end to be zero, this will give a reflection (it is clamped at a not necessarily zero displacement). I think the problem is thinking of this as a wave equation for the displacement of a string (if you are?), I think pressure in a tube with an open (pressure release) end would be more appropriate model to set the boundary condition possibly?.

    CB
    Follow Math Help Forum on Facebook and Google+

  11. #11
    Newbie
    Joined
    Jul 2010
    Posts
    8
    Quote Originally Posted by CaptainBlack View Post
    I think your condition on the supposedly free end is the problem, this looks like the you are constraining the rate of change of U at the end to be zero, this will give a reflection (it is clamped at a not necessarily zero displacement). I think the problem is thinking of this as a wave equation for the displacement of a string (if you are?), I think pressure in a tube with an open (pressure release) end would be more appropriate model to set the boundary condition possibly?.

    CB
    You are probably right :-), but do you have any idea on how can I model it? Actually I'm starting to think that the whole concept of "open end boundary" is wrong for what I need. I've found now an animation on : Open vs Closed pipes (Flutes vs Clarinets)
    and it looks that even if the end is open, there is still some reflection. If this animation is correct, then this mean that my code is fine, but it is not what I need to get. I just want the perturbation to dissapear in "infinity", and not to get any reflection.

    This whole wave analysis problem is a part of my thesis, and analysis of some 3D wave characteristics on reflections on solid boundaries in open space. It works like a charm, until I start to get the reflections from the "walls"(boundaries). I can always make the space bigger, but it is memory consuming. The problem consist in simulating an infinite space in a finite domain.:-)

    How should I proceed?
    Follow Math Help Forum on Facebook and Google+

  12. #12
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by emirs View Post
    You are probably right :-), but do you have any idea on how can I model it? Actually I'm starting to think that the whole concept of "open end boundary" is wrong for what I need. I've found now an animation on : Open vs Closed pipes (Flutes vs Clarinets)
    and it looks that even if the end is open, there is still some reflection. If this animation is correct, then this mean that my code is fine, but it is not what I need to get. I just want the perturbation to dissapear in "infinity", and not to get any reflection.

    This whole wave analysis problem is a part of my thesis, and analysis of some 3D wave characteristics on reflections on solid boundaries in open space. It works like a charm, until I start to get the reflections from the "walls"(boundaries). I can always make the space bigger, but it is memory consuming. The problem consist in simulating an infinite space in a finite domain.:-)

    How should I proceed?
    Since you are using numerical methods to solve this could you not transform the space to a finite line, square, cube,... and then use a zero boundary condition on the boundaries corresponding to the walls at infinity? (this is a top of my head suggestion and the more I think about it the less plausible it seems)

    CB
    Follow Math Help Forum on Facebook and Google+

  13. #13
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    5
    Awards
    2
    In your OP, which derivative were you setting equal to zero at the boundary? Was it U_{t} or U_{x}? Here's a link to something that might be relevant to your computations. I just wonder if you weren't setting the wrong derivative equal to zero. In your discretization, does the first argument correspond to time or space?
    Follow Math Help Forum on Facebook and Google+

  14. #14
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Another probably useless idea: Place the far end of the space grid sufficiently far from the region that you are going to perturb that the signal will not reach it in the time interval you are going to compute over (so if you are going to run the model from 0 to t place the far space boundary out at kct where k is of the order of a few and c the wave propagation speed). Then use a zero boundary condition at the end.

    CB
    Follow Math Help Forum on Facebook and Google+

  15. #15
    Newbie
    Joined
    Jul 2010
    Posts
    8
    Quote Originally Posted by Ackbeet View Post
    In your OP, which derivative were you setting equal to zero at the boundary? Was it U_{t} or U_{x}? Here's a link to something that might be relevant to your computations. I just wonder if you weren't setting the wrong derivative equal to zero. In your discretization, does the first argument correspond to time or space?
    First argument correspond to time (t <-> n), and second to space (x <-> j), so dU_{x}=U(n,j+1)-U(n,j-1), if zero then U(n,j+1)=U(n,j-1) which is used as condition, but doesnt work. I think the free end ("ring" approximation) isn't really what I want to simulate...but to fool a simulator that the string is infinite.
    Follow Math Help Forum on Facebook and Google+

Page 1 of 2 12 LastLast

Similar Math Help Forum Discussions

  1. Closed, open, int, adh and boundary
    Posted in the Differential Geometry Forum
    Replies: 2
    Last Post: March 31st 2011, 05:37 AM
  2. Prove a set is open iff it does not contain its boundary points
    Posted in the Differential Geometry Forum
    Replies: 2
    Last Post: February 24th 2011, 07:16 PM
  3. boundary & initial conditions for wave equation
    Posted in the Differential Equations Forum
    Replies: 2
    Last Post: October 17th 2009, 05:48 PM
  4. Solving a wave equation with inhomogeneous boundary conditions
    Posted in the Differential Equations Forum
    Replies: 2
    Last Post: May 3rd 2009, 03:11 PM
  5. wave equation and Dirichelt boundary condition... HELP!! :(
    Posted in the Advanced Applied Math Forum
    Replies: 1
    Last Post: September 30th 2006, 10:22 PM

Search Tags


/mathhelpforum @mathhelpforum