iconEuler Examples

Surface of an Ellipse

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;

Surface of Ellipsoid

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;

Surface of Ellipsoid

This is about six times as much as the distribution of Thomsens
formula.
>plot3d(x,y,zt-z,points=1,style=".",zoom=4); insimg;

Surface of Ellipsoid

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;

Surface of Ellipsoid

>plot3d(x,y,10*relerror(zt,z),points=1,style=".",zoom=4); insimg;

Surface of Ellipsoid

Examples Homepage