%-------------------------------------------------------------------------------------- % Error curves for the Algebraic Ellipse Perimeter Approximations % and the respective approximations of the elliptic function E(x). %-------------------------------------------------------------------------------------- %-------------------------------------------------------------------------------------- % Globally Best Algebraic Approximations of Ellipse Perimeters %-------------------------------------------------------------------------------------- % This program was written in may 2006 by Stan Sykora while writing % the article on 'Algebraic Approximations of Ellipse Perimeters', % published online at www.ebyte.it/library/docs/math05a/EllipseAlgApprox06.html. % The program is an attachment to the said article and was used to produce Fig.1. %-------------------------------------------------------------------------------------- % The AUTHOR hereby authorizes EVERYBODY to copy, distribute, cut, or modify % the program in whichever form and by whatever means. If you use it and obtain % interesting results, however, you should cite the Web article. % % The AUTOR DISCLAIMS any LIABILITIES which any misguided soul might claim, % due to whatever he/she has done with the program or just thinks to have done. clear;clc;clf;echo off; logscale = 1; % set 0 for linear scale y = logspace(-4,0,10000); y1 = 1+y; [K,E] = ellipke(1-y.*y,1e-15); %Compute the elliptic function E to the highest precision %Actually just 8 digits can be trusted %Cantrell Particularly Fruitful = Rational1, 6313 ppm t1 = (pi-2)/2; % 0.57079632679490 Rational1 = ((1+y.^2)+2*t1*y)./(1+y); plot(y,(Rational1-E)./E,'-k','LineWidth',2,'Color',[0.5 0.3 0]); hold on; err1 = max(abs((Rational1-E)./E)) %Rational2, 557.3 ppm p = [1.22694921875000]; s1 = p(1); t1 = (pi/2)*(1+s1)-1; % 2.49808365277126 Rational2 = ((1+y.^3)+t1*(1+y).*y)./((1+y.^2)+2*s1*y); plot(y,(Rational2-E)./E,'-k','LineWidth',2,'Color',[1 0 0]); hold on; plot(y,-(Rational2-E)./E,'-k','LineWidth',1,'Color',[1 0 0]); hold on; err2 = max(abs((Rational2-E)./E)) %Sykora 2005 = Rational3, 75.96 ppm p = [6.15241658169936 6.16881239339582]; s1 = p(1); t1 = p(2); t2 = (pi/2)*(1+s1)-1-t1; % 4.06617730084445 Rational3 = ((1+y.^4)+t1*y.*(1+y.^2)+2*t2*y.^2)./((1+y.^3)+s1*y.*(1+y)); plot(y,(Rational3-E)./E,'-k','LineWidth',2,'Color',[1 0.7 0]); hold on; plot(y,-(Rational3-E)./E,'-k','LineWidth',1,'Color',[1 0.7 0]); hold on; err3 = max(abs((Rational3-E)./E)) %Rational4, 15.29 ppm p = [13.01750519704827 13.09922140579137 13.02487942169925]; s1 = p(1); s2 = p(2); t1 = p(3); t2 = (pi/2)*(1+s1+s2)-1-t1; % 28.56997512074272 Rational4 = ((1+y.^5)+t1*y.*(1+y.^3)+t2*(y.^2).*(1+y))./((1+y.^4)+s1*y.*(1+y.^2)+2*s2*y.^2); plot(y,(Rational4-E)./E,'-k','LineWidth',2,'Color',[0 0.6 0]); hold on; plot(y,-(Rational4-E)./E,'-k','LineWidth',1,'Color',[0 0.6 0]); hold on; err4 = max(abs((Rational4-E)./E)) %Rational5 = Cantrell 2006, 3.65 ppm p = 100*[0.27297854333670 1.10551872985869 0.27301243680755 1.13302483206429]; s1 = p(1); s2 = p(2); t1 = p(3); t2 = p(4); t3 = (pi/2)*(1+s1+s2)-1-t2-t1; % 76.50091476282086 Rational5 = ((1+y.^6)+t1*y.*(1+y.^4)+t2*(y.^2).*(1+y.^2)+2*t3*(y.^3))./((1+y.^5)+s1*y.*(1+y.^3)+s2*(y.^2).*(1+y)); plot(y,(Rational5-E)./E,'-k','LineWidth',2,'Color',[0 0 0.8]); hold on; plot(y,-(Rational5-E)./E,'-k','LineWidth',1,'Color',[0 0 0.8]); hold on; err5 = max(abs((Rational5-E)./E)) %Rational6, 988 ppb p = 100*[0.51445996674310 3.82256974433855 3.47212007887276 0.51447782789130 3.85327854851892]; s1 = p(1); s2 = p(2); s3 = p(3); t1 = p(4); t2 = p(5); t3 = (pi/2)*(1+s1+s2+s3)-1-t2-t1; % 7.904535392309255e+002 Rational6 = ((1+y.^7)+t1*y.*(1+y.^5)+t2*(y.^2).*(1+y.^3)+t3*(y.^3).*(1+y))./((1+y.^6)+s1*y.*(1+y.^4)+s2*(y.^2).*(1+y.^2)+2*s3*(y.^3)); plot(y,(Rational6-E)./E,'-k','LineWidth',2,'Color',[0 0 0]); hold on; plot(y,-(Rational6-E)./E,'-k','LineWidth',1,'Color',[0 0 0]); hold on; err7 = max(abs((Rational6-E)./E)) %Rational7, 282 ppb p = 1000*[0.09349135794687 1.25936473022183 4.09396201082922 0.09349235523473 1.26273239571330 4.29645229646421]; s1 = p(1); s2 = p(2); s3 = p(3); t1 = p(4); t2 = p(5); t3 = p(6); t4 = (pi/2)*(1+s1+s2+s3)-1-t2-t3-t1; % 2.903735611540449e+003 Rational7 = ((1+y.^8)+t1*y.*(1+y.^6)+t2*(y.^2).*(1+y.^4)+t3*(y.^3).*(1+y.^2)+2*t4*(y.^4))./((1+y.^7)+s1*y.*(1+y.^5)+s2*(y.^2).*(1+y.^3)+s3*(y.^3).*(1+y)); plot(y,(Rational7-E)./E,'-k','LineWidth',2,'Color',[1 0 0]); hold on; plot(y,-(Rational7-E)./E,'-k','LineWidth',1,'Color',[1 0 0]); hold on; err7 = max(abs((Rational7-E)./E)) %-------------------------------------------------------------------------------------- set(gcf,'Color','white'); set(gca,'FontSize',16); set(gca,'LineWidth',2); set(gca,'TickLength',[.03 .03]); set(gca,'yScale','log'); axis([0.0001 1 1e-7 1e-2]); %Below 1e-8 the precision of the program may be insufficient if logscale set(gca,'xScale','log'); %Comment this line for linear x-scale else set(gca,'xTick',[.2,.4,.6,.8,1]); end