distance between two points along a curve
in [Matlab]
|
From: Luca Turchet on 8 May 2010 16:03 Dear Roger, thanks a lot. I tried to understand how that functions proposed by John work, but I was not successiful in applying them for my purposes ;-( I would like to have just an example of the code to use, in order to understand how I can use those functions in my case. For that reason I opted for the solution you proposed. BTW, I implemented it well? It is correct? It is for sure my fault if I have not understood fully how those functions work. I will look at them better now, and I hope to find a solution as soon as possible. Any suggestion or code example is really appreciated ;-) Luca P.S. @John: I forgot to tell you "thanks" for the precious advice. "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hs4eog$dhh$1(a)fred.mathworks.com>... > "Luca Turchet" <tur(a)imi.aau.dk> wrote in message <hs42df$eo0$1(a)fred.mathworks.com>... > > Hi all, > > I think I got it. > > ....... > > I think you should follow John's excellent advice, Luca. It sounds the best to me, since you have a curve y = f(x) with a known formula for its derivative, f'(x). The differential equation you would be solving with 'ode45' is: > > dx/ds = 1/sqrt(1+f'(x).^2) . > > You integrate with respect to arc length s as your independent variable, starting with x as the x-coordinate of point A at s = 0, and you stop when s has reached the distance you wish to travel along the curve to point B. What could be simpler? > > Roger Stafford
From: Roger Stafford on 8 May 2010 17:33 "Luca Turchet" <tur(a)imi.aau.dk> wrote in message <hs4g1p$5jm$1(a)fred.mathworks.com>... > Dear Roger, > thanks a lot. > I tried to understand how that functions proposed by John work, but I was not successiful in applying them for my purposes ;-( > I would like to have just an example of the code to use, in order to understand how I can use those functions in my case. > > For that reason I opted for the solution you proposed. BTW, I implemented it well? > It is correct? > > It is for sure my fault if I have not understood fully how those functions work. I will look at them better now, and I hope to find a solution as soon as possible. > > Any suggestion or code example is really appreciated ;-) > > Luca > > P.S. @John: I forgot to tell you "thanks" for the precious advice. Bruno gave you a good hint as to why your 'find' might come up empty when he said "especially when the function has a large slope". Since you have asked me if your code is correct, perhaps he will forgive me if I spill the beans. The problem is that in the inequalities Z>distance & Z < distance + step_increment you might easily be so unlucky that all Z's are either less than "distance" or greater than "distance+step_increment", since they have larger increments than "step_increment", especially if the slope is still high. Instead of the double inequality, you should be using 'find' with the 'first' option and just your first inequality. Also you need to be sure and select the upper end of the X vector so that it exceeds A_x by at least as much as the quantity "distance" to ensure that a sufficiently large Z will attained. Right now you just have a 15. It should be A_x+distance. One last point. The B_x you have selected gives an arc length beyond the desired distance, and you need to adjust it backwards appropriately. That is, B_x needs to be reduced back towards the previous X value in proportion to the amount back toward the previous Z value of the desired "distance" value. Well, one last, last point. By all means use John's 'ode45' method. It is by far the cleanest of all the ideas that have been put forth here. Don't use mine. It was thought up on the spur of the moment. Roger Stafford
From: John D'Errico on 8 May 2010 21:50 "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hs4lah$756$1(a)fred.mathworks.com>... > Well, one last, last point. By all means use John's 'ode45' method. It is by far the cleanest of all the ideas that have been put forth here. Don't use mine. It was thought up on the spur of the moment. > This is not too difficult to do. Just follow the advice we have offered. For example... Given the function f(x) = sin(x), can I find a point along the curve that lies exactly 25 units in arc length along that curve? I won't do it using a Gaussian, as that is your problem to solve. I will do it several ways though, so you can see the approach. First, I'll use my arc length interpolation routines. I know that if I go out as far as 25 units on the x axis, the arc length along the curve itself must be longer than 25 units. (Any curve must be longer than the straight line distance.) Pick 100 points arbitrarily along that curve. x = linspace(0,25,100); y = sin(x); What is the arc length itself? Use a cubic spline for accuracy here. totallength = arclength(x,y,'spline') totallength = 30.374051851635 It turns out that this curve is a bit longer than the distance along the x axis, but not massively so. Plot the points. plot(x,y,'o') axis equal Now, find that point that is a distance of 25 units along the curve. xy = interparc(25/totallength,x,y,'spline') xy = 20.5883445841352 0.985863637123116 So at (x,y) = (20.5883445841352 , 0.985863637123116) the arc length along the curve should be 25. Of course, this is not exact, since it uses a spline interpolant, then integrates along that spline. We can do this far more accurately by an integration using ode45. The derivative of sin(x) is cos(x). So the differential equation we should solve, parameterized in terms of x, is... dSdx = sqrt(1 + cos(x)^2) Where S(0) = 0. How do we solve it using ode45? fun = @(x,y) sqrt(1 + cos(x)^2); res = ode45(fun,[0,25],0); What is the total arc length of the curve? res.y(end) ans = 30.3786473265496 As you should see, it is quite close to that reported by my own arclength function, which used an integral along a parametric spline curve. But how do we find that point which lies at exactly 25 units along that curve? Easiest is just to start the solver at S(0) = -25. Then I'll just look for a zero crossing. Save this function as an m-file on your search path. % =============================== function [value,isterminal,direction] = ode_events(t,y) % ode event trap, looking for zero crossings of y. value = y; isterminal = ones(size(y)); direction = ones(size(y)); % =============================== Now, call ode45. opts = odeset('events',@ode_events); res = ode45(fun,[0,25],-25,opts); Look at the result. res res = solver: 'ode45' extdata: [1x1 struct] x: [0 2.5 5 7.5 10 12.5 15 17.5 20 20.5870891327006] y: [1x10 double] xe: 20.5870891327006 ye: 7.7715611723761e-16 ie: 1 stats: [1x1 struct] idata: [1x1 struct] The value of res.xe is that point where the arclength to that point is exactly 25 units. Again, it is pretty close to that reported by interparc. The difference is in the approximation of the curve by a spline. Had we sampled the curve at 200 points or more, the error would be smaller. John
From: Luca Turchet on 8 May 2010 21:54 Hi all, I think I have been able to use the John method. Here I post the code. I ask you an opinion. At the end I ask also a MATLAB syntax help: I define a function in the file distance_integral.m (the definition of the files gaussian.m and gaussian_der.m has been given in the previous posts): function [dxds] = distance_integral(s,x) mu = 0;%mu sigma = 2;%sigma dxds = 1/sqrt(1 + gaussian_der(x, mu, sigma).^2); Then I use such function within ODE: %Gaussian parameters mu = 0;%mu sigma = 2;%sigma %Point A coordinates A_x = -0.4; A_y = gaussian(A_x, mu, sigma); distance = 1.2;%wanted distance from A to B Sspan = [0 distance]; % Solve from s=0 to s=distance IC = A_x; % y(s=0) = A_x [S X] = ode45(@distance_integral,Sspan,IC); % Solve ODE B_X = X(length(X)) B_Y = gaussian(B_x, mu, sigma); %PLOT of the gaussian curve starting_point_x = -15; end_point_x = 15; step_increment = 0.001; X = [starting_point_x:step_increment:end_point_x]; Y = gaussian(X, m, s); plot(X,Y) hold on; plot(A_x,A_y,'ro'); plot(B_x,B_y,'k*'); hold off; What do you think? It seems correct to me. Now I have a syntax problem managing the function handles: I would like to pass to the function handled the values of mu and sigma. I tried this but I get an error: [S X] = ode45(@(s,x) distance_integral(s,x,mu,sigma),Sspan,IC) The function in the file distance_integral.m is: function [dxds] = distance_integral(s,x,mu,sigma) dxds = 1/sqrt(1 + gaussian_der(x, mu, sigma).^2); Any suggestion? Please let me know! Cheers "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hs4lah$756$1(a)fred.mathworks.com>... > "Luca Turchet" <tur(a)imi.aau.dk> wrote in message <hs4g1p$5jm$1(a)fred.mathworks.com>... > > Dear Roger, > > thanks a lot. > > I tried to understand how that functions proposed by John work, but I was not successiful in applying them for my purposes ;-( > > I would like to have just an example of the code to use, in order to understand how I can use those functions in my case. > > > > For that reason I opted for the solution you proposed. BTW, I implemented it well? > > It is correct? > > > > It is for sure my fault if I have not understood fully how those functions work. I will look at them better now, and I hope to find a solution as soon as possible. > > > > Any suggestion or code example is really appreciated ;-) > > > > Luca > > > > P.S. @John: I forgot to tell you "thanks" for the precious advice. > > Bruno gave you a good hint as to why your 'find' might come up empty when he said "especially when the function has a large slope". Since you have asked me if your code is correct, perhaps he will forgive me if I spill the beans. The problem is that in the inequalities > > Z>distance & Z < distance + step_increment > > you might easily be so unlucky that all Z's are either less than "distance" or greater than "distance+step_increment", since they have larger increments than "step_increment", especially if the slope is still high. Instead of the double inequality, you should be using 'find' with the 'first' option and just your first inequality. > > Also you need to be sure and select the upper end of the X vector so that it exceeds A_x by at least as much as the quantity "distance" to ensure that a sufficiently large Z will attained. Right now you just have a 15. It should be A_x+distance. > > One last point. The B_x you have selected gives an arc length beyond the desired distance, and you need to adjust it backwards appropriately. That is, B_x needs to be reduced back towards the previous X value in proportion to the amount back toward the previous Z value of the desired "distance" value. > > Well, one last, last point. By all means use John's 'ode45' method. It is by far the cleanest of all the ideas that have been put forth here. Don't use mine. It was thought up on the spur of the moment. > > Roger Stafford
From: Roger Stafford on 8 May 2010 22:53 "John D'Errico" <woodchips(a)rochester.rr.com> wrote in message <hs54cv$3jk$1(a)fred.mathworks.com>... > ........ > dSdx = sqrt(1 + cos(x)^2) > .......... > But how do we find that point which lies at exactly > 25 units along that curve? Easiest is just to start the > solver at S(0) = -25. Then I'll just look for a zero > crossing. Save this function as an m-file on your > search path. > ........ John, why wouldn't it be a lot easier to use 'ode45' to solve the differential equation dxds = 1/sqrt(1+cos(x)^2) as I mentioned earlier? Then you use s as the independent variable in 'ode45' from s = 0 to s = "desired distance", with x starting at the x value of point A. That way you don't have to go to the trouble of searching for a "crossing event" along the way. Roger Stafford
First
|
Prev
|
Next
|
Last
Pages: 1 2 3 4 5 Prev: distance between two points along a curve Next: Matlab Engine Can't Open..... |