Prev: www.jerseysuppliers.com Cheap NFL Jerseys Cheap NHL Jerseys Cheap MLB Jerseys
Next: How to Display Real-time plot(chart), to show process variable value
From: Roland Franzius on 29 May 2010 15:00
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

 | 
Pages: 1
Prev: www.jerseysuppliers.com Cheap NFL Jerseys Cheap NHL Jerseys Cheap MLB Jerseys
Next: How to Display Real-time plot(chart), to show process variable value