Calcul de l'opérateur de courbure


Sur un morceau de surface (mds) régulier x(U) paramétré par la fonction vectorielle x : (u,v) |--> x(u,v) d'un domaine U dans R2 à valeurs dans R3, on considère, en chaque point p, un opérateur A de l'espace tangent Tp qui associe à un vecteur tangent w, la dérivée en p du champ de vecteurs normaux unitaires nu(u,v) (nu = la lettre grecque nu) dans la direction de w.
A est un opérateur auto-adjoint qui, dans la base de Tp formée des vecteurs tangents xu et xv, est représenté par une matrice 2 x 2, disons m.
Cette matrice m est symétrique seulement si cette base est orthonormée. Les invariants de cette matrice, à savoir son déterminant (Det[m]) et sa (demi) trace (1/2 (m[[1,1]] + m[[2,2]])), vont nous permettrent de calculer respectivement la courbure de Gauss K en p et la courbure moyenne H en p (la demi-trace).
La capsule suivante permet le calcul de la matrice m de l'opérateur A en un point quelconque de x(U). Il est conseillé de sauvegarder cette commande pour utilisation future dans un fichier, par exemple, opc.m, dans votre répertoire personnel.

opc[mds_][u_,v_]:= Module[{xu,xv,nu,nor,nuu,nuv,a,b,c,d,syst1,syst2,sol1,sol2},
    xu = D[mds[u,v],u]//Simplify;
    xv = D[mds[u,v],v]//Simplify;
    norme[p_List /; Length[p]==3]:= Sqrt[p . p] //Simplify;
    nu =(1/norme[Cross[xu,xv]]) Cross[xu,xv] // Simplify //PowerExpand ;
    nor[uu_,vv_]:= Evaluate[nu];
    nuu = D[nor[u,v],u]//Simplify;
    nuv = D[nor[u,v],v]//Simplify;
     syst1 = {nuu.xu + a xu.xu + b xv.xu ==0, nuu.xv + a xu.xv + b xv.xv ==0};
     syst2 = {nuv.xu + c xu.xu + d xv.xu ==0, nuv.xv + c xu.xv + d xv.xv ==0};
    sol1 = LinearSolve[{{xu.xu,xv.xu},{xu.xv,xv.xv}},{-nuu .xu,-nuu .xv}];
    sol2 = LinearSolve[{{xu.xu,xv.xu},{xu.xv,xv.xv}},{-nuv .xu,-nuv .xv}];
    m = {sol1,sol2}
    ]
Exemples:
1) Une sphère de rayon r
x[u_,v_]:= {r Cos[u] Sin[v],r Sin[u] Sin[v],r Cos[v]}
opc[x][u,v]
2) Un ellipsoïde
y[u_,v_]:= {2 Cos[u] Sin[v],3 Sin[u] Sin[v],4 Cos[v]}
opc[y][u,v]
Pour cet ellipsoïde, calculons la matrice de l'opérateur de courbure au point (u0,v0) = (.1,.2).
On constate qu'exécuter
opc[y][.1,.2]
ne donne pas le résultat escompté, parce que la commande remplace u par .1 et v par .2 avant de dériver, ce qui engendre un long résultat inutile.
Pour éviter cet obstacle, il suffit d'exécuter
m1 = m /. {u->.1,v->.2}
puisque le résultat de opc[y][u,v] est la matrice m.
On trouve que m1 vaut : {{0.425254, -0.00909543}, {-0.230442, 0.841777}}
Pour afficher cette matrice m1 sous forme de tableau à 2 lignes et 2 colonnes, on exécute
MatrixForm[m1]
Si l'on veut une commande qui accepte comme arguments des valeurs numériques, on peut exécuter la ligne suivante, qui transforme le résultat en une fonction de u et v .
Le nom de cette commande est oc (opérateur de courbure).
On a laissé tomber le p pour éviter une collision avec opc.
oc[mds_][u_,v_]:=( res = opc[mds][s,t];res /. {s->u,t->v})
On conseille d'ajouter la ligne ci-dessus en fin du fichier opc.m pour utilisation ultérieure.
Testons notre nouvelle commande oc sur les exemples précédents :
oc[x][u,v]
oc[x][.1,.2]
oc[y][.1,.1]
La courbure de Gauss K étant le déterminant de la matrice de l'opérateur de courbure, on aura pour cette courbure :
K[mds_][u_,v_]:=Det[ oc[mds][u,v]]
Par exemple la valeur de la courbure de Gauss de l'ellipsoïde au point y(.1,.2) est donnée par
K[y][.1,.2]
on trouve : 0.355874
Il est intéressant de voir le comportement de la courbure de Gauss sur l'ellipsoïde en dessinant le graphe de K.
Il suffit d'exécuter:
Plot3D[K[y][u,v]//Evaluate,{u,0,2Pi},{v,0,Pi}]
Mathematica s'exécute, mais après une liste de messages qu'il est très instructif de lire attentivement. On y découvre qu'il y a un problème lorsque le paramètre v est nul. Et, en effet, notre RPR y n'est pas régulière en tous les points (u,0), comme le montre le calcul du produit vectoriel de xu avec xv :
Les étapes suivantes se passent de commentaires :
norme[v_]:= Sqrt[v .v]
res1 = D[y[u,v],u]
res2 = D[y[u,v],v]
res3 = Cross[res1,res2]
res4 = norme[res3]
res5 = Simplify[PowerExpand[res4]]
InputForm[res5]
grapheradical = Plot3D[35 + 10*Cos[2*u] - 5*Cos[2*(u - v)] -17*Cos[2*v] - 5*Cos[2*(u + v)]//Evaluate,{u,0,2Pi},{v,0,Pi}]
plan0 = Plot3D[0,{u,0,2Pi},{v,0,Pi}]
Options[Plot3D,ViewPoint]
Show[grapheradical,plan0,ViewPoint->{1.3,-2.4,.05}]
Show[grapheradical,plan0,ViewPoint->{1.3,.5,.05}]
Plot3D[K[y][u,v]//Evaluate,{u,0,2Pi},{v,.1,3.14}]

Retour à ma page personnelle



©jmt : création le 001101 - dernière modification le 011221