Post by Roger StaffordWell, 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