From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: upper bound of integration (numerical evaluation)
Date: Wed, 25 Mar 2009 19:06:02 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 21
Message-ID: <gqdvaq$66$>
References: <gq8jfd$c3j$> <> <gqalvi$ge8$> <gqb4cp$rkj$> <gqdqa1$tb$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: 1238007962 198 (25 Mar 2009 19:06:02 GMT)
NNTP-Posting-Date: Wed, 25 Mar 2009 19:06:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:527706

"Dimitar Dimitrov" <> wrote in message <gqdqa1$tb$>...
> 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