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,
(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:
and leads to:
Unfortunately this condition creates a partial reflection from the end which is supposed to be free.
Please help
Regards
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.
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:
(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.
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
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
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?
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
In your OP, which derivative were you setting equal to zero at the boundary? Was it or ? 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?
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
First argument correspond to time (t <-> n), and second to space (x <-> j), so , 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.