by R. Grothmann The following function computes the surface of an ellipse with semiaxis a, b, and c. We use elliptical integrals to evaluate the surface.
>function map SE (a,b,c) ...
v=sort([a,b,c]); a=v[3]; b=v[2]; c=v[1]; if c~=0 then return b*2*pi; endif; s=sqrt(a^2-c^2); if s~=0 then return 4*pi*a^2; endif; t=acos(c/a); m=(b^2-c^2)/(b*sin(t))^2; return 2*pi*(c^2+b*s*elle(t,m)+b*c^2/s*ellf(t,m)); endfunction
We wish to fit a function to these data. To this end, we evaluate SE(1,x,y) on a data set with spacing 1/n, and x>=y. First we generate this data set using the Euler matrix language. First we generate a matrix of x and y data, then we flatten out thes matrices to a row vector.
>n=100; t=linspace(0,1,n); >x=redim(dup(t,n+1),1,(n+1)^2); >y=redim(dup(t',n+1),1,(n+1)^2);
Now we use only the data we want to have.
>i=nonzeros(x>=y); x=x[i]; y=y[i];
Then we compute the surface on these points. This function takes about 10 seconds on my computer. The elliptical integrals are implemented in Euler in the Euler language using Carlson's iteration.
>z=SE(1,x,y);
Now we can plot the data.
>plot3d(x,y,z/50,points=1,style=".",angle=-60°); insimg;
Knud Thomsen proposed the following approximation. He supposes the parameter p=1.6075. We adjust p here, so that the formula is exact for c=0, and for a=b=c.
>function SEThomsen (a,b,c,p=log(3)/log(2)) ...
return 4*pi*(((a*b)^p+(b*c)^p+(a*c)^p)/3)^(1/p) endfunction
>zt=SEThomsen(1,x,y);
The maximal relative error is about 1.4%.
>totalmax(abs(relerror(zt,z)))
0.0139929973544
Let us try to minimize the parameter p for our test grid.
>function SEtest(p) ...
global x,y,z; zt=SEThomsen(1,x,y,p); return totalmax(abs(relerror(zt,z))); endfunction
The minimal p is close to 1.6, which is the proposed value of Thomsen.
>fmin("SEtest",1,2)
1.60729709574
Next, we try to fit a quadratic function to the data. These functions are of the form bx+cy+dxy+ex^2+fy^2. To find the coefficiens, we have to generate a matrix M, containing the values of these basis functions at our points.
>M=(x_y_x*y_x^2_y^2)'; b=z';
We fit the minimal relative error, so we minimze |M/b.a-1|.
>afit=fit(M,b)
5.96124797529 2.77365954996 0.555431568826 0.178478223732 3.19578868983
The maximal error is more than 0.10.
>totalmax(M.afit-b)
0.104337631173
The relative error is more than 11%.
>totalmax(abs(relerror(M.afit,b)))
0.111031287968
So we try a best approximation fit. First we set up a function, which computes the maximal relative error for specific coefficients. Note that the Nelder-Mead method needs a row vector.
>function minhelp(a) := totalmax(abs(M.a'-b));
This Nelder-Mead method takes about one second on my computer.
>amin=neldermin("minhelp",afit')'
6.06715135423 2.69752900417 0.586679981676 0.117535805136 3.1867265446
But the maximal absolute error did not change much.
>minhelp(amin')
0.0992689312478
In the following plot, we can see the distribution of the absolute error of this formula.
>plot3d(x,y,(M.amin-z')',>points,style=".",zoom=4); insimg;
This is about six times as much as the distribution of Thomsens formula.
>plot3d(x,y,zt-z,points=1,style=".",zoom=4); insimg;
Next, we try to minimize the relative error.
>function minhelprel(a) := totalmax(abs(relerror(M.a',b)));
Now, the Nelder-Mead method takes a little bit longer.
>arelmin=neldermin("minhelprel",amin')'
6.05936287466 3.47275247797 -1.91629733585 0.40301203644 4.6732247235
The maximal relativ error is now 3.4%
>totalmax(abs(relerror(M.arelmin,b)))
0.0349810329329
>plot3d(x,y,10*relerror(M.arelmin,z')',points=1,style=".",zoom=4); insimg;
>plot3d(x,y,10*relerror(zt,z),points=1,style=".",zoom=4); insimg;