Surfaces de niveaux avec scilab

Tracer une représentation graphique des niveaux d'une fonction quelconque n'est pas un problème anodin les quelques fois où je tombe sur cette question je ne peux m'en sortir que sur des cas particulier où l'on peut séparer les variables. Pour une fonction à deux variables $f:{\mathbb R}^2\longmapsto {\mathbb R}$, suffisamment régulière, les niveaux $S=\{(x,y)\vert f(x,y)=C\}$ sont des courbes du plan,  Scilab possède des fonctions prédéfinies contour et contourf pour tracer ces courbes de niveaux. Mais ces fonctions ne peuvent pas traiter le cas de fonctions à trois variables $f:{\mathbb R}^3\longrightarrow{\mathbb R}$, où les niveaux $S=\{(x,y,z)\vert f(x,y,z)=C\}$ sont des surfaces dans l'espace ... heureusement j'ai trouvé un code intéressant sur Scilab File Exchange définissant une fonction contour3d capable de tracer des surfaces de niveaux juste à partir de la fonction  f et de la valeur  C, petit aperçu en vidéo :

La séquence d'appel  de la fonction   est contour3d(x,y,z,f,nf) 

• x,y,z  sont des vecteurs de taille nx,ny,nz qui contiennent le découpage  une discrétisation des valeurs des trois axes de coordonnées et définissent donc une grille  dans le cube où l'on cherche les niveaux de la fonction
• f est une hypermatrice de taille nx x ny x nz contenant les valeurs de la fonction sur la grille définie par x,y,z
• nf est un vecteur contenant les valeurs des niveaux que l'on cherche à tracer

pour détecter les niveaux des fonctions  il faut commencer par  les  évaluer sur un maillage cubique. Pour les fonctions à une ou deux variables on dispose de la fonction scilab feval  mais celle-ci ne s'applique pas aux fonctions à trois variables! Pour se faciliter la tâche le plus simple est de  créer une  fonction équivalente à feval qui à partir d'une discrétisation de trois intervalles x,y,z à n1,n2,n3 valeurs, va construire une hypermatrice de taille n1xn2xn3  contenant toute les valeurs de la fonction f sur ce maillage :

function t=feval3d(x,y,z,f)
    n1=length(x);n2=length(y);;n3=length(z);// nombre de points 
    t=zeros(n1,n2,n3)//hypermatrice  pour contenir les valeurs de f(x,y,z)
    for i=1:n1
        for j=1:n2
            for k=1:n3
                t(i,j,k)=f(x(i),y(j),z(k));
            end
        end
    end
endfunction

On peut maintenant expliquer  les exemples  de la vidéo, ils sont basé sur les deux fonctions $f_1(x,y,z)=x^2+y^2+z^2$  et  $f_2(x,y,z)=x^2+y^2-z^2$ :

function t=f1(x,y,z)
    t=x.^2+y.^2+z.^2
endfunction

function t=f2(x,y,z)
    t=x.^2+y.^2-z.^2
endfunction

• pour la fonction f1 les niveaux pour C=1,2,3,4 sont des sphères emboîtées :
  
  
  

  // discrétisation du domaine
  x=[-2.5:0.2:2.5];y=x;z=x; 
  t=feval3d(x,y,z,f1);
  //affichage de 4 niveaux  
  clf
  contour3d(x,y,z,t,[1:4]);
  

• pour la fonction f2  les niveaux sont des hyperboloïdes , à une nappe  si C>0 à deux nappes si C<0 et un cône si C=0 :
  
  
  

  // discrétisation du domaine
  x=[-2:0.2:2];y=x;z=x; 
  t=feval3d(x,y,z,f2);
  //affichage de 3 niveaux 
  clf
  contour3d(x,y,z,t,[-1 0 2]); 
  

• mais on peut aussi tracer sur une même figure les niveaux de deux fonctions différentes, dans ce cas il faut modifier les couleurs des surfaces manuellement pour qu'on puisse les différencier :
  
  
  

  x=[-2.5:0.2:2.5];y=x;z=x; 
  t1=feval3d(x,y,z,f1);
  t2=feval3d(x,y,z,f2);
  clf
  contour3d(x,y,z,t1,1); 
  E=gce();E.color_mode=20; //modif couleur surface 
  contour3d(x,y,z,t2,-1);
  E=gce();E.color_mode=10;// modif couleur surface 
  

le calcul des  niveaux peut nécessiter pas mal de mémoire pour les calculs il peut être utile d'étendre la mémoire disponible pour scilab avec stacksize('max') avant de réaliser ces exemples. Pour ceux que ça intéresse voici le code légèrement modifié (entre les lignes 62 et 76) que j'ai utilisé :

function contour3d(x,y,z,f,nf,varargin)
// Draws 3D contours, iso-surfaces, of function f on
// 3D lattice x,y,z.
// DESCRIPTION contour3d  draws level surfaces of a
// 3D function f(x,y,z) on a 3D plot.  The values of
// f(x,y,z) are to be in the 3D hypermatrix f at the
// lattice grid points defined by 1D or 3D matrices 
// x, y and z.
// 
// INPUT  If no arguments are supplied at all, then
// contour3d draws a demonstration of nested spheres.
// x,y,z: 1D or 3D matrices of coordinates of grid
// points in space at which the function f is known.
// f: 3D matrix of function values. Set any missing
// f-values to %nan.
// nf: if nf is simply a number of surfaces to draw,
// then a colorbar is also drawn.  However, if nf is
// a vector of surface values, then no colorbar. You
// can skip one or more colors in the current
// colormap by specifying %nan at the corresponding
// place in the vector nf of f values---useful for
// superimposing iso-surface plots.
// 
// Any extra arguments to contour3d are passed on to
// plot3d to empower you to set view, legend, and so on.
// 
// OUTPUT fvals: vector of set contour values.
// 
// Uses 40.nx.ny.nz memory or thereabouts.
// Tony Roberts, 26 Aug 2006; modified Mar 2010.
// Philippe Roux  modified 1 July 2013
[nlhs,nrhs]=argn()
if nrhs==0,
  nx=3;ny=4;nz=5;
  x=linspace(0,1,nx);
  y=linspace(0,1,ny);
  z=linspace(0,1,nz);
  // make row/column/plane of hypermat
  xx=hypermat([nx,1,1],x);
  yy=hypermat([1,ny,1],y);
  zz=hypermat([1,1,nz],z);
  xx=xx(:,ones(1,ny),ones(1,nz));
  yy=yy(ones(1,nx),:,ones(1,nz));
  zz=zz(ones(1,nx),ones(1,ny),:);
  // compute function and pass on to draw
  f=sqrt((xx-1).^2+yy.^2+zz.^2)-0.5;
  nf=6;
  varargin=list();
  clf()
end

// first get & check sizes of 3D arrays
if length(size(f))~=3, error("contour3d: f must be 3D hypermatrix"), end
if min(size(f))<2 data-blogger-escaped-contour3d:="" data-blogger-escaped-error="" data-blogger-escaped-need="">=2 points in each dimension"), end
[nx,ny,nz]=size(f);
if length(x(:))==nx, x=hypermat([nx,1,1],x(:)); x=x(:,ones(1,ny),ones(1,nz)); end
if length(y(:))==ny, y=hypermat([1,ny,1],y(:)); y=y(ones(1,nx),:,ones(1,nz)); end
if length(z(:))==nz, z=hypermat([1,1,nz],z(:)); z=z(ones(1,nx),ones(1,ny),:); end
if size(x)~=[nx,ny,nz],  error("contour3d: x 3D matrix does not match f"), end
if size(y)~=[nx,ny,nz],  error("contour3d: y 3D matrix does not match f"), end
if size(z)~=[nx,ny,nz],  error("contour3d: z 3D matrix does not match f"), end
// find iso-surface values to include zero
F=gcf();F.color_map=hotcolormap(32);
if (length(nf)==1)  // only one iso-surface
  makefvals=%F;// no colorbar  -->  //minf=nf-1;maxf=nf+1;  df=1
  fvals=nf
  color_flag=0// use color_mode  to change color
else // multiple iso-surface (nf is a vector)
  makefvals=%T;
  fvals=gsort(nf,'g','i')
  F.color_map=hotcolormap(length(nf)+1);
  maxf=max(nf); minf=min(nf);
  df=(maxf-minf)/(length(nf)-1+%eps);
  maxf=maxf+df; minf=minf;//
  color_flag=2// 
end
// make a huge list of tetrahedrons
// use these to cater for 8 different orientations of the 5 tetras
i=2*ceil((1:nx-1)/2); id=(-1).^(1:nx-1);
j=2*ceil((1:ny-1)/2); jd=(-1).^(1:ny-1);
k=2*ceil((1:nz-1)/2); kd=(-1).^(1:nz-1);
l=(nx-1)*(ny-1)*(nz-1); // number of 3D bricks
ntet=5*l;  // number of tetras
// define tetras by their lattice point labels
p=hypermat([nx,ny,nz],1:nx*ny*nz); // index of points
ps=zeros(ntet,4);
ps(:,1)=matrix([p(i,j,k);p(i+id,j,k+kd);p(i,j+jd,k+kd);p(i+id,j+jd,k);p(i+id,j,k)],ntet,1);
ps(:,2)=matrix([p(i+id,j,k);p(i,j,k+kd);p(i,j,k+kd);p(i,j+jd,k);p(i,j+jd,k)],ntet,1);
ps(:,3)=matrix([p(i,j+jd,k);p(i+id,j,k);p(i,j+jd,k);p(i+id,j,k);p(i,j,k+kd)],ntet,1);
ps(:,4)=matrix([p(i,j,k+kd);p(i+id,j+jd,k+kd);p(i+id,j+jd,k+kd);p(i+id,j+jd,k+kd);p(i+id,j+jd,k+kd)],ntet,1);
fs=[f(ps(:,1)) f(ps(:,2)) f(ps(:,3)) f(ps(:,4))];
// sort each tetra corner into decreasing order of f
[fs,ll]=gsort(fs,"c");
ps(:)=ps(sub2ind([ntet,4],(1:ntet)'*ones(1,4),ll));
// only use tetras that have all vertices defined
ok=find(~isnan(sum(fs,"c")));

// build facets for each isosurface value
xf=[]; yf=[]; zf=[]; cf=[];
for col=1:length(fvals)
  fval=fvals(col)
  i=sum(fs(ok,:)>fval,"c"); // how many values>contour
  // do all the 1-3 tetras
  j=ok(find(i==1)); 
  if ~isempty(j) then k=ones(1,4); l=[2 2:4];
  rats=(fval-fs(j,l))./(fs(j,k)-fs(j,l));
  xf=[xf; rats.*matrix(x(ps(j,k)),-1,4)+(1-rats).*matrix(x(ps(j,l)),-1,4)];
  yf=[yf; rats.*matrix(y(ps(j,k)),-1,4)+(1-rats).*matrix(y(ps(j,l)),-1,4)];
  zf=[zf; rats.*matrix(z(ps(j,k)),-1,4)+(1-rats).*matrix(z(ps(j,l)),-1,4)];
  cf=[cf col(ones(j))];
  end
  // do all the 2-2 tetras
  j=ok(find(i==2)); 
  if ~isempty(j) then k=[1 1 2 2]; l=[3 4 4 3];
  rats=(fval-fs(j,l))./(fs(j,k)-fs(j,l));
  xf=[xf; rats.*matrix(x(ps(j,k)),-1,4)+(1-rats).*matrix(x(ps(j,l)),-1,4)];
  yf=[yf; rats.*matrix(y(ps(j,k)),-1,4)+(1-rats).*matrix(y(ps(j,l)),-1,4)];
  zf=[zf; rats.*matrix(z(ps(j,k)),-1,4)+(1-rats).*matrix(z(ps(j,l)),-1,4)];
  cf=[cf col(ones(j))];
  end
  // do all the 3-1 tetras
  j=ok(find(i==3)); 
  if ~isempty(j) then k=[1 1:3]; l=4*ones(1,4);
  rats=(fval-fs(j,l))./(fs(j,k)-fs(j,l));
  xf=[xf; rats.*matrix(x(ps(j,k)),-1,4)+(1-rats).*matrix(x(ps(j,l)),-1,4)];
  yf=[yf; rats.*matrix(y(ps(j,k)),-1,4)+(1-rats).*matrix(y(ps(j,l)),-1,4)];
  zf=[zf; rats.*matrix(z(ps(j,k)),-1,4)+(1-rats).*matrix(z(ps(j,l)),-1,4)];
  cf=[cf col(ones(j))];
  end
end
// draw a colorbar only if nf is a vector of requested iso-surfaces
if makefvals, colorbar(minf,maxf), end
// draw facets according to current colormap
plot3d(xf',yf',list(zf',cf'),varargin(:))
fog=gce();
fog.hiddencolor=-1; // make front and back same color
fog.color_flag=color_flag;
endfunction