Calcul de la courbure géodésique d'une courbe sur une surface

La courbure géodésique kappag d'une courbe donnée par une représentation naturelle sur une surface S en un de ses points P est la norme de la projection de son vecteur de courbure en P sur le plan tangent à S en P. Si x(u(s),v(s)) est cette courbe, alors le calcul de sa courbure géodésique revient au calcul du produit mixte de x', x'' et nu, et fait intervenir, via les équations de Gauss, les symboles de Christoffel et les dérivées u',v' et u'', v''.
La capsule ci-dessous va grandement faciliter le calcul de cette courbure dans Mathematica. On conseille de la sauvegarder pour utilisation ultérieure dans un fichier (je suggère kappag.m) de son répertoire personnel.
Contenu du fichier kappag.m :
(*-- calcul de la courbure géodésique d'une courbe sur ---*)
(*-- un morceau de surface mds ---*)
kappag[mds_][u_,v_][s_]:= Module[{x, a, b, der1, der2, EE, FF, GG, disc, EEu, EEv, FFu, FFv, GGu, GGv, c11, c12, c22, kg, p, q},
    x[a_,b_]:= mds[a,b];
    der1 = D[x[a,b],a];
    der2 = D[x[a,b],b];
    EE= der1 . der1;
    FF= der1 . der2;
    GG= der2 . der2;
    disc = EE GG -FF^2;
    EEu = D[EE,a];
    EEv = D[EE,b];
    FFu = D[FF,a];
    FFv = D[FF,b];
    GGu = D[GG,a];
    GGv = D[GG,b];
    c11 =1/( 2 disc) {EEv*FF - 2*FF*FFu + EEu*GG,-EE*EEv - EEu*FF + 2*EE*FFu};
    c12 =1/( 2 disc) {EEv*GG - FF*GGu,-EEv*FF + EE*GGu};
    c22 =1/( 2 disc){2*FFv*GG - GG*GGu - FF*GGv,-2*FF*FFv + FF*GGu + EE*GGv};
    m = {{c11[[1]],c11[[2]]},{c12[[1]],c12[[2]]},{c22[[1]],c22[[2]]}} /. Thread[Rule[{a,b},{u,v}]];
    disc = disc /.Thread[Rule[{a,b},{u,v}]];
    kg = Sqrt[disc] ( m[[1,2]] p^3 + (2 m[[2,2]] - m[[1,1]])p^2 q + (m[[3,2]] - 2 m[[2,1]]) p q^2 - m[[3,1]] q^3 + p q' - q p');
    subst = Thread[Rule[{u,v,p,q},{u[s],v[s],u'[s],v'[s]}]];
    kapg = kg /. subst
    ]
Il faudra encore simplifier au besoin …
Il est bon de tester sa commande sur quelques exemples :
N'oublions pas qu'avec la commande kappag, nous calculons la courbure géodésique d'une courbe quelconque (paramétrée par l'abcisse curviligne) sur un morceau de surface quelconque.
Il faut donc s'attendre à ce que le résultat soit long, et il est recommandé, comme d'habitude, de nommer les objets qu'on calcule pour pouvoir y référer facilement.

x[u_,v_] := {u,v,u+v}
res1 = kappag[x][u,v][s]

y[u_,v_] := {u,v,u^2 + v^2}
res2 = kappag[y][u,v][s]
res3 = Simplify[res2]

Pour l'ellipsoïde général, donné par :

ell[u_,v_]:= {a Cos[u] Sin[v],b Sin[u] Sin[v],c Cos[v]}
res4 = kappag[ell][u,v][s]
On peut chercher à simplifier ce long résultat, mais sans grand succès, à cause de la présence des constantes a, b et c, non nécessairement égales.
Au lieu de faire le calcul de la courbure géodésique d'une courbe sur un ellipsoïde général, spécialisons le résultat obtenu pour le cas particulier d'une sphère de rayon R en remplaçant, dans le résultat 4 (res4), les constantes a, b, c par R :
res5 = res4 /. Thread[Rule[{a,b,c},{R,R,R}]]
Ce faisant, nous allons rejoindre nos préoccupations des premières leçons du cours où nous avons considéré des grands cercles sur la terre comme sphère de rayon R, en particulier des grands cercles passant par Montréal.
Après quelques simplifications, on considère la courbure géodésique de l'équateur (elle est nulle, donc l'équateur est une géodésique), celle d'un cercle parallèle passant par Montréal (non nulle!), celle du méridien passant par Montréal (nulle, car grand cercle!) :
courbgeodsphere = Simplify[res5]
courbgeodsphere = PowerExpand[%]
courburegeodesiqueequateur = courbgeodsphere /. Thread[ Rule[{u[s], v[s], u'[s], v'[s], u''[s], v''[s]},{s/R,Pi/2,1/R,0,0,0}]]
courburegeodesiqueparallele = courbgeodsphere /. Thread[ Rule[{u[s], v[s], u'[s], v'[s], u''[s], v''[s]},{Sqrt[2] s/R, Pi/4,Sqrt[2]/R ,0,0,0}]]
courburegeodesiquemeridien = courbgeodsphere /. Thread[ Rule[{u[s], v[s], u'[s], v'[s], u''[s], v''[s]},{0,s/R,0, 1/R,0,0}]]

Retour à ma page personnelle



©jmt : création le 001129 - dernière modification le 011121