|
|
Re: ellipse points regularly spaced, or inverse of the incomplete elliptic integral of the second (not first) kind
Posted:
Jun 1, 2010 3:45 PM
|
|
"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
|
|
