Thread Subject:

upper bound of integration (numerical evaluation)

Hi,
I have a curve parametrized using a parameter "T".
The relation between "T" and the length of the curve "L" is expressed through an integral of some function "F", let us denote it with L = quad(F,x_i,x_f)
where x_i and x_f are lower and upper bounds for the integral.

If we have a given value for L, I would like to find such x_f that produces it.

The problem can be easily solved using bisection, however, I was wondering whether there is a standard Matlab function that can be used.

Thanks,
Dimitar

"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gq8gvc$cps$1@fred.mathworks.com>...
> Hi,
> I have a curve parametrized using a parameter "T".
> The relation between "T" and the length of the curve "L" is expressed through an integral of some function "F", let us denote it with L = quad(F,x_i,x_f)
> where x_i and x_f are lower and upper bounds for the integral.
>
> If we have a given value for L, I would like to find such x_f that produces it.
>
> The problem can be easily solved using bisection, however, I was wondering whether there is a standard Matlab function that can be used.
>

There is no standard function for your explicit
problem, but fzero is a rootfinder which will
be a better choice than writing your own
bisection code.

John

>
> There is no standard function for your explicit
> problem, but fzero is a rootfinder which will
> be a better choice than writing your own
> bisection code.
>
> John

Hi John,
Thank you for the post!
I am not sure how to use fzero in order to solve my problem.
Could you elaborate a bit...

Regards,
Dimitar

> >
> > There is no standard function for your explicit
> > problem, but fzero is a rootfinder which will
> > be a better choice than writing your own
> > bisection code.
> >
> > John
>
> Hi John,
> Thank you for the post!
> I am not sure how to use fzero in order to solve my
> problem.
> Could you elaborate a bit...
>
> Regards,
> Dimitar

Determine the zero of the function
f(x) = L - int_{x_i}^{x} F(t) dt
by using fzero.

Best wishes
Torsten.

> Determine the zero of the function
> f(x) = L - int_{x_i}^{x} F(t) dt
> by using fzero.
>
> Best wishes
> Torsten.

of course ...
Thanks,
Dimitar

"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gqalvi$ge8$1@fred.mathworks.com>...
> > Determine the zero of the function
> > f(x) = L - int_{x_i}^{x} F(t) dt
> > by using fzero.
> >
> > Best wishes
> > Torsten.
>
> of course ...
> Thanks,
> Dimitar

  A difficulty with using 'fzero' in a straightforward manner for your kind of problem is that with each call on the objective function you would define, it will commence an integration again all the way back from the original lower integration limit. This means the same territory is being integrated over and over again as 'fzero' makes its numerous repeated calls. Even though 'fzero' may be very efficient as to the total number of calls it makes, this strikes me as an inherently inefficient procedure.

  You might consider defining a special objective function that stores a cumulative integral value and at each call simply algebraically adds to that value by starting its next lower integration limit from the previous upper limit value, even if it is necessary to integrate backwards. This way the unnecessary repetition in continually starting back at the original lower limit could be much reduced.

  The function 'fzero' is presumably designed to minimize the total number of calls it makes on a user-defined objective function. However in your particular problem it is not so much the total number of calls made but the total absolute length of integration path covered that should be minimized. For this reason you might reconsider your earlier idea of some kind of sectioning procedure combined with the above cumulative technique but designed to minimize the total integration path length likely to be covered.

Roger Stafford

Hi Roger,
Thanks for the comment.

Of course, thats how I implemented it, otherwise (as you said) the solver is integrating over the same region again and again.

Even though fzero implements a nice algorithm, I do not think that it is the best choice for my problem. First, because I can give not only a guess x0 that is close to the solution, but an interval [a b] around the solution. Hence, I don't need the "Phase I" (of determining [a b]) performed by fzero. And second, the function in this case permits using directly a "secant" algorithm which converges faster than a bisection or standard linear interpolation.

Of course, other options exist, but I expected that for such a standard problem there would already be a worked solution in Matlab. I say standard, because for example if one wants to relate the length of a given curve to some other parameterization "t" of the same curve (which is what I do), leads to this problem.

Bye the way, my curve is defined as a piecewise cubic polynomial, and coming to think of it, probably there is a completely alternative way to approach this problem. Do you have any suggestions?

Cheers,
Dimitar

"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gqdqa1$tb$1@fred.mathworks.com>...
> Hi Roger,
> Thanks for the comment.
>
> Of course, thats how I implemented it, otherwise (as you said) the solver is integrating over the same region again and again.
>
> Even though fzero implements a nice algorithm, I do not think that it is the best choice for my problem. First, because I can give not only a guess x0 that is close to the solution, but an interval [a b] around the solution. Hence, I don't need the "Phase I" (of determining [a b]) performed by fzero. And second, the function in this case permits using directly a "secant" algorithm which converges faster than a bisection or standard linear interpolation.
>
> Of course, other options exist, but I expected that for such a standard problem there would already be a worked solution in Matlab. I say standard, because for example if one wants to relate the length of a given curve to some other parameterization "t" of the same curve (which is what I do), leads to this problem.
>
> Bye the way, my curve is defined as a piecewise cubic polynomial, and coming to think of it, probably there is a completely alternative way to approach this problem. Do you have any suggestions?
>
> Cheers,
> Dimitar

  I neglected to mention it before but there is a matlab function called 'cumtrapz' which could prove useful to you here. (Locating the upper limit from the output of 'cumtrapz' is simply a matter of searching its vector output for the point where the value 'L' is crossed.) By itself it may not be accurate enough for your purposes but used in conjunction with 'quad' it might be competitive with the "secant" method you mentioned. Used once over a sufficiently long interval to include the desired upper limit and with a reasonably fine resolution, 'cumtrapz' could give you a fairly good approximation of the right interval to try and 'quad' could then make this more accurate within that interval. Then another 'cumtrapz' confined to this interval followed by 'quad' for greater accuracy would pin things down even closer. I am not sure how many such repetitions are necessary but it
shouldn't be very many.

  As for piecewise cubic polynomials, within any cubic segment, arclength is some kind of elliptic integral (I forget the details) and its inverse would therefore be an elliptic function. Matlab has functions which directly compute elliptic functions, called 'ellipj' and 'ellipke', and which could presumably solve your problem in one step with no iterations once you have identified the proper segment.

Roger Stafford

Hi Roger,
It has been a while since I used elliptic integrals.
Do you have a reference where I can find a way to represent the standard arch length integral

int( dx/dt^2 + dy/dt^2 + dz/dt^2 )_{a}^{b}

with an elliptic integral (if I remember correctly it should be of the second kind).

Today I read some articles regarding my problem, and there is no mention about such possibility (only alternative numerical methods are discussed). My "guess" is that there should be a good reason for that (though I have no idea what it is)! It would be really educational if someone can clarify this point.

By the way, I am quite satisfied with the performance of the "secant" method. It converges in average using less that five function evaluations, even though my interval [a b] is [x_current t(end)] (which can clearly be improved :)).

Cheers,
Dimitar

"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gqefm1$lj9$1@fred.mathworks.com>...
> It has been a while since I used elliptic integrals.
> Do you have a reference where I can find a way to represent the standard arch length integral
>
> int( dx/dt^2 + dy/dt^2 + dz/dt^2 )_{a}^{b}
> .......

  I am in the same boat as you are, Dimitar. It has been a long time since I did anything with either elliptic integrals or elliptic functions. I did well to remember that elliptic functions were inverses of certain elliptic integrals. I don't have any good references except the obvious textbooks on the subject and my textbooks are all very old.

  By the way, in the integral you wrote you need to take the square root of that integrand to get a valid arc-length:

 int(sqrt(dx/dt^2+dy/dt^2+..))dt, ...

Roger Stafford

> > int( dx/dt^2 + dy/dt^2 + dz/dt^2 )_{a}^{b}

Well, it was a bit late in the evening ...

Thanks to all for the help. I will stick to my numerical solution (unless someone intervenes).

Regards,
Dimitar