Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
ellipse points regularly spaced, or inverse of the incomplete

Subject: ellipse points regularly spaced, or inverse of the incomplete

From: Felipe G. Nievinski

Date: 26 May, 2010 17:30:01

Message: 1 of 9

Are you aware of an inverse of the incomplete elliptic integral of the
second (not first) kind? I'm familiar with the Jacobi amplitude as the
inverse of the incomplete elliptic integral of the first kind, but
can't seem to find in the literature an inverse for the second kind.
Alternatively, it'd help to know how to obtain the first kind given
the second one (both incomplete).

I found this statement at <http://dlmf.nist.gov/19.25#v>:
    "[elliptic] integrals of the second kind, ... are not invertible
in terms of single-valued functions"
though I find the treatment given, in terms of the so-called elliptic
symmetric integrals, a little hard to digest.

I'm trying to discretize the perimeter of an ellipse at regularly
spaced points.Trivial discretizations end up causing an undesirable
high/low density of points at maximum/minimum distance from the
center. (A brute force approach would be to discretize at a
sufficiently high density everywhere then discard excessive points
where necessary.)

I am given the ellipse semi-major and semi-minor axes, a and b,
respectively; eccentricity is obtained as e=sqrt(1-(b/a)^2). To define
regularly spaced points, I compute the ellipse perimeter or
circumference length C and, given a step length s, define intermediary
arc lengths as len = 0:s:C. Now the problem is to obtain the x,y
coordinates of the ellipse at each such arc length values. (The
circumference is given by C = 4 * a * E(k), where E(k) is the complete
elliptical integral of the second kind, and k=e^2).

A brief review of the background theory. The ellipse arc length can be
calculated as:
  len = a * E(phi,k),
in terms of the incomplete elliptical integral of the _second_ kind
  v = E(phi,k).
I'm not aware of an explicit formulation for its inverse,
  phi = Einv(v,k).
In contrast, for the incomplete elliptical integral of the _first_
kind,
  u = F(phi,k)
its inverse is the so-called Jacobi's amplitude function:
  phi = am(u,k),

Sometimes the arc length is expressed as:
  len = a * Eps(u,k),
in terms of the so-called Jacobi's epsilon function:
  Eps(u,k) = E(am(u,k),k).
I don't see the benefit of this formulation, unless we can relate the
output of the first and second elliptic integrals, u and v. And if we
do, Jacobi's amplitude function usually involves the so-called Jacobi
elliptic functions:
  sn(u) = sin(phi),
  cn(u) = cos(phi),
which can be related directly to the point position by
  x = a * sn(u,k),
  y = b * cn(u,k).

Putting it all together, we could proceed in one of the following
ways:
  v = len / a
  invert for phi = Einv(v,k) (how?)
  x = a * cos(phi)
  y = b * sin(phi)
or
  v = len / a
  relate u to v (how?)
  get sn(u,k) and cn(u,k) given u via routine ellipj.m
  x = a * sn(u,k)
  y = b * cn(u,k)

I'm inclined to obtain phi = Einv(v,k) via a semi-brute force
approach, i.e., sampling densely v = E(phi,k) then keeping only the
phi values yielding more or less regularly spaced v, e.g.:
  % (using <http://elliptic.googlecode.com/svn/trunk/elliptic12.m>)
  s = 0.1;
  a = 1;
  b = 1/2;
  esq = 1 - (b/a)^2;
  phi = linspace(0, 2*pi)';
  [F, E] = elliptic12 (phi, esq);
  len = a * E;
  dlen = [0; diff(len)];
  figure, plotyy(phi,len, phi,dlen)
  [ignore, temp] = ellipke (esq);
  C = 4 * a * temp;
    C - len(phi==2*pi) % DEBUG
  len2 = (0:s:C)'; len2(end+1) = C;
  phi2 = interp1(len, phi, len2, 'nearest');
  dlen2 = [0; diff(len2)];
  figure, plotyy(phi2,len2, phi2,dlen2)
  figure, hold on, plot(phi,len,'.-k'), plot(phi2,len2,'ok')
  figure, hold on, plot(phi,dlen,'.-k'), plot(phi2,dlen2,'ok')
  %x = a * cos(phi);
  %y = b * sin(phi);
  varphi = atan2(a*sin(phi), b*cos(phi));
  x = a * cos(varphi);
  y = b * sin(varphi);
  ind = ismember(phi, phi2);
  figure, hold on, plot(x, y, '.-k'), plot(x(ind), y(ind), 'ok'), axis
equal
  figure, hold on, plot(x(ind), y(ind), 'o-k'), axis equal

Though if I'm going this way I might just as well avoid integrals
altogether and do simply a pure brute force approach:
  s = 0.1;
  a = 1;
  b = 1/2;
  theta = linspace(0, 2*pi)';
  x = a * cos(theta);
  y = b * sin(theta);
  figure, plot(x, y, '.-k'), axis equal
  dx = [0; diff(x)];
  dy = [0; diff(y)];
  dlen = sqrt(dx.^2 + dy.^2);
  len = cumsum(dlen);
  figure, plotyy(theta,len, theta,dlen), title('len')
  C = len(phi==2*pi);
  len2 = (0:s:C)'; len2(end+1) = C;
  theta2 = interp1(len, theta, len2, 'nearest');
  dlen2 = [0; diff(len2)];
  figure, plotyy(theta2,len2, theta2,dlen2), title('E')
  figure, hold on, plot(theta,dlen,'.-k'), plot(theta2,dlen2,'ok')
  ind = ismember(theta, theta2);
  figure, hold on, plot(x, y, '.-k'), plot(x(ind), y(ind), 'ok'), axis
equal
  figure, hold on, plot(x(ind), y(ind), 'o-k'), axis equal

References:
<http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html>
<http://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html>
<http://mathworld.wolfram.com/JacobiAmplitude.html>
<http://mathworld.wolfram.com/JacobiEllipticFunctions.html>
<http://mathworld.wolfram.com/Ellipse.html>
<http://en.wikipedia.org/wiki/Ellipse>
<http://en.wikipedia.org/wiki/Jacobi's_elliptic_functions>
<http://en.wikipedia.org/wiki/Elliptic_integral>
<http://reference.wolfram.com/mathematica/tutorial/
EllipticIntegralsAndEllipticFunctions.html>
<http://www.mhtl.uwaterloo.ca/courses/me755/web_chap3.pdf>
<http://dlmf.nist.gov/19>
<http://dlmf.nist.gov/22>

Subject: ellipse points regularly spaced, or inverse of the incomplete

From: Bruno Luong

Date: 26 May, 2010 19:19:22

Message: 2 of 9

"Felipe G. Nievinski" <fgnievinski@gmail.com> wrote in message <htjlqp$e3q$1@news.acm.uiuc.edu>...

>
> I found this statement at <http://dlmf.nist.gov/19.25#v>:
> "[elliptic] integrals of the second kind, ... are not invertible
> in terms of single-valued functions"
> though I find the treatment given, in terms of the so-called elliptic
> symmetric integrals, a little hard to digest.
>

Are they talking about invertible in C or in R? To illustrate, y = x^3 is invertible in R but not in C. This is surely an "incomplete" statement. LOL

Your post is too long to read, and did not even spot any question. Is it just a discussion?

Bruno

Subject: ellipse points regularly spaced, or inverse of the incomplete

From: John D'Errico

Date: 26 May, 2010 19:41:04

Message: 3 of 9

"Felipe G. Nievinski" <fgnievinski@gmail.com> wrote in message <htjlqp$e3q$1@news.acm.uiuc.edu>...
> Are you aware of an inverse of the incomplete elliptic integral of the
> second (not first) kind? I'm familiar with the Jacobi amplitude as the
> inverse of the incomplete elliptic integral of the first kind, but
> can't seem to find in the literature an inverse for the second kind.
> Alternatively, it'd help to know how to obtain the first kind given
> the second one (both incomplete).
>
> I found this statement at <http://dlmf.nist.gov/19.25#v>:
> "[elliptic] integrals of the second kind, ... are not invertible
> in terms of single-valued functions"
> though I find the treatment given, in terms of the so-called elliptic
> symmetric integrals, a little hard to digest.
>
> I'm trying to discretize the perimeter of an ellipse at regularly
> spaced points.Trivial discretizations end up causing an undesirable
> high/low density of points at maximum/minimum distance from the
> center. (A brute force approach would be to discretize at a
> sufficiently high density everywhere then discard excessive points
> where necessary.)

Better, and quite simple, is just to generate a set of
arbitrarily spaced points on the ellipse, and then use
my interparc function, to interpolate points equally
spaced on the ellipse.

http://www.mathworks.com/matlabcentral/fileexchange/27096-interparc

This will be reasonably fast, and quite accurate.

HTH,
John

Subject: ellipse points regularly spaced, or inverse of the incomplete

From: John D'Errico

Date: 26 May, 2010 20:00:24

Message: 4 of 9

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <htjtgg$gvj$1@fred.mathworks.com>...
> "Felipe G. Nievinski" <fgnievinski@gmail.com> wrote in message <htjlqp$e3q$1@news.acm.uiuc.edu>...
> > Are you aware of an inverse of the incomplete elliptic integral of the
> > second (not first) kind? I'm familiar with the Jacobi amplitude as the
> > inverse of the incomplete elliptic integral of the first kind, but
> > can't seem to find in the literature an inverse for the second kind.
> > Alternatively, it'd help to know how to obtain the first kind given
> > the second one (both incomplete).
> >
> > I found this statement at <http://dlmf.nist.gov/19.25#v>:
> > "[elliptic] integrals of the second kind, ... are not invertible
> > in terms of single-valued functions"
> > though I find the treatment given, in terms of the so-called elliptic
> > symmetric integrals, a little hard to digest.
> >
> > I'm trying to discretize the perimeter of an ellipse at regularly
> > spaced points.Trivial discretizations end up causing an undesirable
> > high/low density of points at maximum/minimum distance from the
> > center. (A brute force approach would be to discretize at a
> > sufficiently high density everywhere then discard excessive points
> > where necessary.)
>
> Better, and quite simple, is just to generate a set of
> arbitrarily spaced points on the ellipse, and then use
> my interparc function, to interpolate points equally
> spaced on the ellipse.
>
> http://www.mathworks.com/matlabcentral/fileexchange/27096-interparc
>
> This will be reasonably fast, and quite accurate.
>

Alternatively, you could also just use the same scheme I
use in that tool, here applying an odesolver (ode45) to
the ellipse, parameterized by angle around the center.

HTH,
John

Subject: ellipse points regularly spaced, or inverse of the incomplete

From: Roger Stafford

Date: 27 May, 2010 02:56:04

Message: 5 of 9

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <htjuko$1ug$1@fred.mathworks.com>...
> Alternatively, you could also just use the same scheme I
> use in that tool, here applying an odesolver (ode45) to
> the ellipse, parameterized by angle around the center.
>
> HTH,
> John

  Applying John's idea, you could use arclength, s, as the independent variable in ode45 and solve the simultaneous equations

 dx/ds = -a^2*y/sqrt(b^4*x^2+a^4*y^2)
 dy/ds = +b^2*x/sqrt(b^4*x^2+a^4*y^2)

to travel counterclockwise around the ellipse. You need only go a quarter of the way around. You set the s steps returned and stopping point in accordance with the perimeter values you already know.

Roger Stafford

Subject: ellipse points regularly spaced, or inverse of the incomplete elliptic integral of the second (not first) kind

From: David W. Cantrell

Date: 27 May, 2010 19:00:01

Message: 6 of 9

"Felipe G. Nievinski" <fgnievinski@gmail.com> wrote:
> Are you aware of an inverse of the incomplete elliptic integral of the
> second (not first) kind? I'm familiar with the Jacobi amplitude as the
> inverse of the incomplete elliptic integral of the first kind, but can't
> seem to find in the literature an inverse for the second kind.
> Alternatively, it'd help to know how to obtain the first kind given the
> second one (both incomplete).

I recommend that you look at the sci.math.num-analysis thread "Dividing
an ellipse into equal parts" and particularly at my post therein of
5 Mar 2006:

<http://groups.google.com/group/sci.math.num-analysis/msg/1c3cd4eb79b8e86e>

It might give what you need.

[By the way, the second link mentioned there, to the sci.math.research
thread "Inverting elliptic integrals", seems to be broken. However, that
thread can now be found at the Math Forum:
<http://mathforum.org/kb/message.jspa?messageID=1672061>.]

Regards,
David W. Cantrell

> I found this statement at <http://dlmf.nist.gov/19.25#v>:
> "[elliptic] integrals of the second kind, ... are not invertible
> in terms of single-valued functions"
> though I find the treatment given, in terms of the so-called
> elliptic symmetric integrals, a little hard to digest.
>
> I'm trying to discretize the perimeter of an ellipse at regularly
> spaced points.Trivial discretizations end up causing an undesirable
> high/low density of points at maximum/minimum distance from the
> center. (A brute force approach would be to discretize at a
> sufficiently high density everywhere then discard excessive points
> where necessary.)

[snip]

Subject: ellipse points regularly spaced, or inverse of the incomplete

From: Roland Franzius

Date: 29 May, 2010 19:00:02

Message: 7 of 9

Felipe G. Nievinski schrieb:
> Are you aware of an inverse of the incomplete elliptic integral of the
> second (not first) kind? I'm familiar with the Jacobi amplitude as the
> inverse of the incomplete elliptic integral of the first kind, but
> can't seem to find in the literature an inverse for the second kind.
> Alternatively, it'd help to know how to obtain the first kind given
> the second one (both incomplete).
>
> I found this statement at <http://dlmf.nist.gov/19.25#v>:
> "[elliptic] integrals of the second kind, ... are not invertible
> in terms of single-valued functions"
> though I find the treatment given, in terms of the so-called elliptic
> symmetric integrals, a little hard to digest.
>
> I'm trying to discretize the perimeter of an ellipse at regularly
> spaced points.Trivial discretizations end up causing an undesirable
> high/low density of points at maximum/minimum distance from the
> center. (A brute force approach would be to discretize at a
> sufficiently high density everywhere then discard excessive points
> where necessary.)
>
> I am given the ellipse semi-major and semi-minor axes, a and b,
> respectively; eccentricity is obtained as e=sqrt(1-(b/a)^2). To define
> regularly spaced points, I compute the ellipse perimeter or
> circumference length C and, given a step length s, define intermediary
> arc lengths as len = 0:s:C. Now the problem is to obtain the x,y
> coordinates of the ellipse at each such arc length values. (The
> circumference is given by C = 4 * a * E(k), where E(k) is the complete
> elliptical integral of the second kind, and k=e^2).
>
> A brief review of the background theory. The ellipse arc length can be
> calculated as:
> len = a * E(phi,k),
> in terms of the incomplete elliptical integral of the _second_ kind
> v = E(phi,k).
> I'm not aware of an explicit formulation for its inverse,
> phi = Einv(v,k).
> In contrast, for the incomplete elliptical integral of the _first_
> kind,
> u = F(phi,k)
> its inverse is the so-called Jacobi's amplitude function:
> phi = am(u,k),
>
> Sometimes the arc length is expressed as:
> len = a * Eps(u,k),
> in terms of the so-called Jacobi's epsilon function:
> Eps(u,k) = E(am(u,k),k).
> I don't see the benefit of this formulation, unless we can relate the
> output of the first and second elliptic integrals, u and v. And if we
> do, Jacobi's amplitude function usually involves the so-called Jacobi
> elliptic functions:
> sn(u) = sin(phi),
> cn(u) = cos(phi),
> which can be related directly to the point position by
> x = a * sn(u,k),
> y = b * cn(u,k).
>
> Putting it all together, we could proceed in one of the following
> ways:
> v = len / a
> invert for phi = Einv(v,k) (how?)
> x = a * cos(phi)
> y = b * sin(phi)
> or
> v = len / a
> relate u to v (how?)
> get sn(u,k) and cn(u,k) given u via routine ellipj.m
> x = a * sn(u,k)
> y = b * cn(u,k)
>
> I'm inclined to obtain phi = Einv(v,k) via a semi-brute force
> approach, i.e., sampling densely v = E(phi,k) then keeping only the
> phi values yielding more or less regularly spaced v, e.g.:
> % (using <http://elliptic.googlecode.com/svn/trunk/elliptic12.m>)
> s = 0.1;
> a = 1;
> b = 1/2;
> esq = 1 - (b/a)^2;
> phi = linspace(0, 2*pi)';
> [F, E] = elliptic12 (phi, esq);
> len = a * E;
> dlen = [0; diff(len)];
> figure, plotyy(phi,len, phi,dlen)
> [ignore, temp] = ellipke (esq);
> C = 4 * a * temp;
> C - len(phi==2*pi) % DEBUG
> len2 = (0:s:C)'; len2(end+1) = C;
> phi2 = interp1(len, phi, len2, 'nearest');
> dlen2 = [0; diff(len2)];
> figure, plotyy(phi2,len2, phi2,dlen2)
> figure, hold on, plot(phi,len,'.-k'), plot(phi2,len2,'ok')
> figure, hold on, plot(phi,dlen,'.-k'), plot(phi2,dlen2,'ok')
> %x = a * cos(phi);
> %y = b * sin(phi);
> varphi = atan2(a*sin(phi), b*cos(phi));
> x = a * cos(varphi);
> y = b * sin(varphi);
> ind = ismember(phi, phi2);
> figure, hold on, plot(x, y, '.-k'), plot(x(ind), y(ind), 'ok'), axis
> equal
> figure, hold on, plot(x(ind), y(ind), 'o-k'), axis equal
>
> Though if I'm going this way I might just as well avoid integrals
> altogether and do simply a pure brute force approach:
> s = 0.1;
> a = 1;
> b = 1/2;
> theta = linspace(0, 2*pi)';
> x = a * cos(theta);
> y = b * sin(theta);
> figure, plot(x, y, '.-k'), axis equal
> dx = [0; diff(x)];
> dy = [0; diff(y)];
> dlen = sqrt(dx.^2 + dy.^2);
> len = cumsum(dlen);
> figure, plotyy(theta,len, theta,dlen), title('len')
> C = len(phi==2*pi);
> len2 = (0:s:C)'; len2(end+1) = C;
> theta2 = interp1(len, theta, len2, 'nearest');
> dlen2 = [0; diff(len2)];
> figure, plotyy(theta2,len2, theta2,dlen2), title('E')
> figure, hold on, plot(theta,dlen,'.-k'), plot(theta2,dlen2,'ok')
> ind = ismember(theta, theta2);
> figure, hold on, plot(x, y, '.-k'), plot(x(ind), y(ind), 'ok'), axis
> equal
> figure, hold on, plot(x(ind), y(ind), 'o-k'), axis equal
>
> References:
> <http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html>
> <http://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html>
> <http://mathworld.wolfram.com/JacobiAmplitude.html>
> <http://mathworld.wolfram.com/JacobiEllipticFunctions.html>
> <http://mathworld.wolfram.com/Ellipse.html>
> <http://en.wikipedia.org/wiki/Ellipse>
> <http://en.wikipedia.org/wiki/Jacobi's_elliptic_functions>
> <http://en.wikipedia.org/wiki/Elliptic_integral>
> <http://reference.wolfram.com/mathematica/tutorial/
> EllipticIntegralsAndEllipticFunctions.html>
> <http://www.mhtl.uwaterloo.ca/courses/me755/web_chap3.pdf>
> <http://dlmf.nist.gov/19>
> <http://dlmf.nist.gov/22>
>


EllipticE(u,k) is a linear function in u plus a logarithmic derivative
of the elliptic theta_4-series, which starts with 1 -2 cos (... u).

So a fourier approximation with one two terms should be sufficient in
most cases, if the ellipse is not extremely excentric.

See eg Gradshteyn/Rhyzik 8.194

E(u,k)= u(1- Th''(0)/Th(0)) + Th'(u)/Th(u)

and the formulas defining Th, K , K' and so on.

Be aware that the modulus k of Gradshteyn is k^2 in Abramowitz/Mathematica.

Once you have determined the period for perodic parts cos(n pi u/K), the
fourier transforms should give nice approximations with few terms.

--

Roland Franzius

Subject: ellipse points regularly spaced, or inverse of the incomplete

From: Felipe G. Nievinski

Date: 1 Jun, 2010 17:30:02

Message: 8 of 9

Thanks a lot for all your suggestions! Now I have to digest them. ;)
-Felipe.
PS: it came to my attention that not all the posts are showing up at
the Google usenet mirror,
  <http://groups.google.com/group/sci.math.num-analysis/browse_thread/
thread/e4f1b49c83189281/7fdd52d226c31096>;
you can see all the msgs (currently seven) at
  <http://www.mathworks.com/matlabcentral/newsreader/view_thread/
283132>;
thanks to moiseev.igor for point it out.

Subject: ellipse points regularly spaced, or inverse of the incomplete elliptic integral of the second (not first) kind

From: Bruno Luong

Date: 1 Jun, 2010 19:45:21

Message: 9 of 9

"David W. Cantrell" <DWCantrell@sigmaxi.net> wrote in message <htmffh$jnc$1@news.acm.uiuc.edu>...
>
> [By the way, the second link mentioned there, to the sci.math.research
> thread "Inverting elliptic integrals", seems to be broken. However, that
> thread can now be found at the Math Forum:
> <http://mathforum.org/kb/message.jspa?messageID=1672061>.]
>
> Regards,
> David W. Cantrell

Incidentally, I bump to the same problem today, and I try David's method. Unless if I make an error somewhere, the power series seem to have poor accuracy for m ~ 1 and z >> 0.

I hope I make a mistake somewhere:

function E = elliptic12_powerseries(z, m)
% Compute EllipticE by power series in z
% Provided by David Cantrell
% http://mathforum.org/kb/message.jspa?messageID=1672061

% E = z + ...
% (m*z^3)/6 + ...
% (1/120)*(-4*m + 13*m^2)*z^5 + ...
% ((16*m - 284*m^2 + 493*m^3)*z^7)/5040 + ...
% ((-64*m + 4944*m^2 - 31224*m^3 + 37369*m^4)*z^9)/362880 + ...
% ((256*m - 81088*m^2 + 1406832*m^3 - 5165224*m^4 +
% 4732249*m^5)*z^11)/39916800;

P = {1;
     (1/6)*[1 0];
     (1/120)*[13 -4 0];
     (1/5040)*[493 -284 16 0];
     (1/362880)*[37369 -31224 4944 -64 0];
     (1/39916800)*[4732249 -5165224 1406832 -81088 256 0] };
P = cellfun(@(p) polyval(p,m), P);
E = z.*polyval(P(end:-1:1),z.^2);

end % elliptic12_powerseries

%%
m = 0.9;
theta=linspace(0,pi/2);
Epowerseries = elliptic12_powerseries(theta, m);
% http://www.mathworks.com/matlabcentral/fileexchange/8805
[trash Eanalytic] = elliptic12(theta, m);

plot(theta, Eanalytic, 'r', theta, Epowerseries, 'b')
title(sprintf('Ellipictic integral 2nd kind, m = %g', m));
xlabel('z (rad)')
legend('Analytic', 'power series','Location','NorthWest');

Any comment?

Bruno

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us