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 SCILAB, 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