Bonjour, j'ai fais plusieurs codes sur les "CARACTÉRISATION DE SURFACE PAR LA MÉTHODE DU PONT CAPILLAIRE" en MATLAB,

j'aimerais savoir si quelqu'un pourrait traduire cela en SCILAB, svp !

La partie la plus dure est celle de la résolution d'équation différentielle de Laplace (du 2nd ordre). Savez vous résoudre celle ci sur SICLAB, sans passer par des solvers svp ?

Pour plus d'information vous pouvez m'envoyer un mail : janot92@hotmail.fr

Voici ci dessous les codes :


Fonction Dérivé de Y par rapport à z :

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
function [ v ] = DYdz( z,Y ) %v : paramètre de sortie ; z et y : paramètres d'entrée 
global lc; %lc = kappa 
v=[Y(2) 
  (1+Y(2)^2)/Y(1)+(1+Y(2)^2)^(3/2)*z/(lc^2)]; 
end
Fonction Point capillaire (en 2D) :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
function [z,Y] = ptcap2D (h,theta_deg) 
global lc; 
warning off; 
lc=2; 
R=100; 
z_lentille=linspace(0,10,101); 
r_lentille=sqrt(2*R.*(z_lentille-h)-((z_lentille-h).^2)); 
zpas=1; 
zini=h+0.01; 
for iii=1:1:4 
   zpas=zpas/(10); 
   for z0=zini:zpas:R/10 
       r0=sqrt(2*R.*(z0-h)-((z0-h)^2)); 
       alpha=asin(r0/R); 
       theta=theta_deg/180*pi; 
       rp0=cot(alpha+theta); 
       %solveur 
       [z,Y]=ode45('DYdz',[z0 0],[r0 rp0]); 
       taille=size(z); 
       zfinal=z(taille(1)); 
       taille2=size(Y); 
       y2final=Y(2,taille2(2)); 
       if (zfinal~=0 || abs(y2final)>100) 
           zini=z0-zpas; 
           [z,Y]=ode45('DYdz',[zini 0],[r0 rp0]); 
           break 
       end 
   end 
end 
%tracés (Attention : il faut mettre ces tracés en commentaire pour obtenir les autres courbes) 
hold on; 
 
plot(r_lentille,z_lentille); 
plot(-r_lentille,z_lentille); 
plot(Y(:,1),z); 
plot(-Y(:,1),z); 
axis equal; 
axis([-60 60 0 10]) 
end 
Fonction Pont capillaire (en 3D) : 
function [X,Y,Zlentille] = ptcap3D ( ) 
R=100; 
h=0; 
for zi=0:1:100 
      r_lentille=sqrt(2.*R.*(zi/10-h)-(zi/10-h).^2); 
   for the=0:100 
      X(zi+1,the+1)=r_lentille.*cos(the*3.6*pi/180); 
      Y(zi+1,the+1)=r_lentille.*sin(the*3.6*pi/180); 
   end 
end 
for xi=1:101 
   for yi=1:101 
      Zlentille(xi,yi)=((2*R+2*h)... 
         -sqrt(((-2*R-2*h)^2)-4.*(2*R*h+h^2+X(xi).^2+Y(yi).^2)))./2; 
   end 
end 
[Zpt,Ytmp]=ptcap2D(h,0);            % on se sert de ptcap2D 
tailleY=size(Ytmp); 
for zi=1:tailleY(1,1) 
   r_pt=Ytmp(zi,1); 
   for the=1:tailleY(1,1) 
      Xpt(zi,the)=r_pt.*cos((the-1)/tailleY(1,1)*360*pi/180); 
      Ypt(zi,the)=r_pt.*sin((the-1)/tailleY(1,1)*360*pi/180); 
   end 
end 
for xi=1:tailleY(1,1) 
   for yi=1:tailleY(1,1) 
      Zpt_tab(xi,yi)=Zpt(xi); 
   end 
end 
%tracés du graphe 3D 
figure( ) 
hold on 
mesh(X,Y,Zlent); % la function mess équivaut à un plot, mais en 3D 
mesh(Xpt,Ypt,Zpt_tab); 
colormap cool; %colormap est une fonction permettant d'avoir plusieurs coloris du graphe 
axis equal; 
 
hold off; 
end
Fonction variation de Thêta (pour obtenir la courbe A0=f(), pour h=0) :

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
function variation_theta () 
h = 0; 
kkk = 0; 
for theta_deg=0:1:178 
   kkk = kkk +1; 
   [z,Y] =ptcap2D(h,theta_deg); 
   tab_theta(kkk,1)=theta_deg; % on stock les valeurs de theta_deg dans un tableau 
   Aire(kkk,1)=pi*Y(1,1)^2; 
   %pi*Y(1,1)^2 
end 
hold on 
plot(tab_theta,Aire,,'+r'); 
hold off 
end
Fonction variation de h (pour obtenir la courbe A=f(h) pour différents  fixé) :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
function variation_h() 
num_theta=0 ; 
for theta_deg=0:20:180 
   %faire varier h avec un thêta fixe qui se trouve dans argument grâce à un for 
% ici on varie thêta entre 0° et 180° 
num_theta= num_theta+1 ; 
hpas=0.5; 
   %Aire=zeros(100,1); 
   %tab_h=zeros(100,1); 
   jjj=0; 
   for jjj=0:1:1 
       kkk=0; 
       for h=0:hpas:100 
          kkk=kkk+1; 
          [z,Y] =ptcap2D(h,theta_deg); 
          %h 
          tab_h(kkk, num_theta)=h; 
          Aire(kkk, num_theta)=pi*Y(1,1)^2; 
          %pi*Y(1,1)^2 
          taille=size(z); 
          zfinal=z(taille(1)); 
          if(zfinal~=0) 
              if(kkk<=30) 
 
             hpas=hpas/10; 
             break; 
          else 
             break; 
          end 
       end 
        %utiliser ptcap2D pour avoir dernière valeur du tableau z qui correspond au pont capillaire 
       %si cette valeur est différente de 0 on arrête avec un break 
       %On doit prendre la première valeur de Y pour calculer l'aire de contact 
       %pi*r2=pi*Y2 => aire de contact 
       %tracer a(h) 
       if(kkk>30) 
          break; 
       end 
    end 
  end 
  hold on 
  plot(tab_h,Aire,'+r'); 
  hold off 
end 
end