http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330
MATLAB Central Newsreader  upper bound of integration (numerical evaluation)
Feed for thread: upper bound of integration (numerical evaluation)
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Mon, 23 Mar 2009 17:30:20 +0000
upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637006
Dimitar Dimitrov
Hi, <br>
I have a curve parametrized using a parameter "T". <br>
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)<br>
where x_i and x_f are lower and upper bounds for the integral.<br>
<br>
If we have a given value for L, I would like to find such x_f that produces it. <br>
<br>
The problem can be easily solved using bisection, however, I was wondering whether there is a standard Matlab function that can be used. <br>
<br>
Thanks,<br>
Dimitar

Mon, 23 Mar 2009 17:46:01 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637016
John D'Errico
"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gq8gvc$cps$1@fred.mathworks.com>...<br>
> Hi, <br>
> I have a curve parametrized using a parameter "T". <br>
> 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)<br>
> where x_i and x_f are lower and upper bounds for the integral.<br>
> <br>
> If we have a given value for L, I would like to find such x_f that produces it. <br>
> <br>
> The problem can be easily solved using bisection, however, I was wondering whether there is a standard Matlab function that can be used. <br>
> <br>
<br>
There is no standard function for your explicit<br>
problem, but fzero is a rootfinder which will<br>
be a better choice than writing your own<br>
bisection code.<br>
<br>
John

Mon, 23 Mar 2009 18:13:01 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637023
Dimitar Dimitrov
> <br>
> There is no standard function for your explicit<br>
> problem, but fzero is a rootfinder which will<br>
> be a better choice than writing your own<br>
> bisection code.<br>
> <br>
> John<br>
<br>
Hi John,<br>
Thank you for the post! <br>
I am not sure how to use fzero in order to solve my problem.<br>
Could you elaborate a bit...<br>
<br>
Regards,<br>
Dimitar

Tue, 24 Mar 2009 07:24:58 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637157
Torsten Hennig
> > <br>
> > There is no standard function for your explicit<br>
> > problem, but fzero is a rootfinder which will<br>
> > be a better choice than writing your own<br>
> > bisection code.<br>
> > <br>
> > John<br>
> <br>
> Hi John,<br>
> Thank you for the post! <br>
> I am not sure how to use fzero in order to solve my<br>
> problem.<br>
> Could you elaborate a bit...<br>
> <br>
> Regards,<br>
> Dimitar<br>
<br>
Determine the zero of the function<br>
f(x) = L  int_{x_i}^{x} F(t) dt<br>
by using fzero.<br>
<br>
Best wishes<br>
Torsten.

Tue, 24 Mar 2009 13:08:02 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637226
Dimitar Dimitrov
> Determine the zero of the function<br>
> f(x) = L  int_{x_i}^{x} F(t) dt<br>
> by using fzero.<br>
> <br>
> Best wishes<br>
> Torsten.<br>
<br>
of course ...<br>
Thanks, <br>
Dimitar

Tue, 24 Mar 2009 17:14:02 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637317
Roger Stafford
"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gqalvi$ge8$1@fred.mathworks.com>...<br>
> > Determine the zero of the function<br>
> > f(x) = L  int_{x_i}^{x} F(t) dt<br>
> > by using fzero.<br>
> > <br>
> > Best wishes<br>
> > Torsten.<br>
> <br>
> of course ...<br>
> Thanks, <br>
> Dimitar<br>
<br>
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.<br>
<br>
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.<br>
<br>
The function 'fzero' is presumably designed to minimize the total number of calls it makes on a userdefined 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.<br>
<br>
Roger Stafford

Wed, 25 Mar 2009 17:40:17 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637673
Dimitar Dimitrov
Hi Roger, <br>
Thanks for the comment.<br>
<br>
Of course, thats how I implemented it, otherwise (as you said) the solver is integrating over the same region again and again.<br>
<br>
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. <br>
<br>
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.<br>
<br>
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?<br>
<br>
Cheers,<br>
Dimitar

Wed, 25 Mar 2009 19:06:02 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637697
Roger Stafford
"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gqdqa1$tb$1@fred.mathworks.com>...<br>
> Hi Roger, <br>
> Thanks for the comment.<br>
> <br>
> Of course, thats how I implemented it, otherwise (as you said) the solver is integrating over the same region again and again.<br>
> <br>
> 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. <br>
> <br>
> 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.<br>
> <br>
> 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?<br>
> <br>
> Cheers,<br>
> Dimitar<br>
<br>
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 <br>
shouldn't be very many. <br>
<br>
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.<br>
<br>
Roger Stafford

Wed, 25 Mar 2009 23:45:05 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637775
Dimitar Dimitrov
Hi Roger,<br>
It has been a while since I used elliptic integrals. <br>
Do you have a reference where I can find a way to represent the standard arch length integral <br>
<br>
int( dx/dt^2 + dy/dt^2 + dz/dt^2 )_{a}^{b}<br>
<br>
with an elliptic integral (if I remember correctly it should be of the second kind).<br>
<br>
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.<br>
<br>
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 :)). <br>
<br>
Cheers,<br>
Dimitar

Thu, 26 Mar 2009 00:50:03 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637788
Roger Stafford
"Dimitar Dimitrov" <mail_mitko@example.com> wrote in message <gqefm1$lj9$1@fred.mathworks.com>...<br>
> It has been a while since I used elliptic integrals. <br>
> Do you have a reference where I can find a way to represent the standard arch length integral <br>
> <br>
> int( dx/dt^2 + dy/dt^2 + dz/dt^2 )_{a}^{b}<br>
> .......<br>
<br>
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.<br>
<br>
By the way, in the integral you wrote you need to take the square root of that integrand to get a valid arclength:<br>
<br>
int(sqrt(dx/dt^2+dy/dt^2+..))dt, ...<br>
<br>
Roger Stafford

Thu, 26 Mar 2009 06:38:03 +0000
Re: upper bound of integration (numerical evaluation)
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247330#637826
Dimitar Dimitrov
> > int( dx/dt^2 + dy/dt^2 + dz/dt^2 )_{a}^{b}<br>
<br>
Well, it was a bit late in the evening ... <br>
<br>
Thanks to all for the help. I will stick to my numerical solution (unless someone intervenes). <br>
<br>
Regards,<br>
Dimitar