Géodésique sur une sphère, passant par un point donné, dans une direction donnée.

Une courbe sur une surface est une géodésique si, en chacun de ses points, sa courbure géodésique est nulle. On sait que le vecteur de courbure géodésique est la projection, sur le plan tangent à la surface, du vecteur de courbure de la courbe en chacun de ses points. Il s'en suit qu'un grand cercle de la sphère (un cercle dans un plan passant par le centre de la sphère) est une géodésique.
Dans cette page, nous allons voir comment déterminer, sur l'exemple simple de la sphère, la géodésique passant par un point donné et dans une direction donnée.
Cet exemple va contenir toutes les étapes nécessaires et va permettre d'appliquer la méthode utilisée à des surfaces plus compliquées.
Il s'agira de résoudre NUMERIQUEMENT une équation différentielle ordinaire NON LINEAIRE du deuxième ordre, traduisant la condition que la courbure géodésique de la courbe cherchee est identiquement nulle.
Soit donc une représentation paramétrique réguliere de la sphere de rayon 1, centrée à l'origine :
x[u_,v_]:= {Cos[u] Sin[v],Sin[u] Sin[v],Cos[v]} (* -Pi < u < Pi, 0 < v < Pi *)
Dessinons-la :
sphere = ParametricPlot3D[x[u,v]//Evaluate,{u,-Pi,Pi},{v,0,Pi}]
Choisissons le point (0,Pi/4) dans l'espace des paramètres et une direction t dans cet espace, p.ex. parallèle à la bissectrice du premier quadrant lorsque t = 1.
Chargeons la commande kappag pour calculer la courbure géodésique :
<<kappag.m
Nommons les résultats de nos calculs pour référence ultérieure :
res1 = kappag[x][u,v][s]
res2 = Simplify[PowerExpand[res1]]
res3 = PowerExpand[res2]
Dans ces calculs, s est l'abscisse curviligne de la courbe cherchée.
Choisisons la paramétrisation suivante pour la courbe dans l'espace des paramètres (u,v) dont l'image devra être une géodésique : u(s) = s et v(s) = une fonction v de s à trouver, telle que la courbe x(s,v(s)) soit une géodésique.
Avec les choix faits plus haut, on aura donc à procéder aux substitutions suivantes dans res3 : u[s] -> s, u'[s] -> 1, u''[s] -> 0
res4 = res3 /. {u[s]->s,u'[s]->1,u''[s]->0}
L'expression obtenue, posée égale à zéro, est l'équation différentielle qui va nous permettre de trouver la géodésique cherchée satisfaisant aux conditions initiales. On a :
equadiff = res4 == 0
La solution en fonction de la direction t est fournie par une intégration numérique grâce à la commande suivante :
sol[t_]:= NDSolve[{equadiff,v[0]==Pi/4,v'[0]==t},v,{s,-1,1}]
On rappelle que DSolve tente d'intégrer des équations ou systèmes d'équations différentielles symboliquement et que, par contre, NDSolve utilise des algorithmes numériques et donc ne tolère pas d'arguments symboliques.
La commande sol[t] ci-dessus n'a qu'un but : rendre l'écriture du problème plus pratique. Il suffit, par exemple, d'exécuter sol[1] si l'on veut la géodésique dans la direction v'[0] = 1, ou sol[.3] dans la direction de pente .3.
Avant d'exécuter sol[1], sachons que Mathematica va nous donner le résultat sous la forme d'une règle de substitution de v par une fonction numérique "InterpolatingFunction" , définie sur le domaine {s,-1,1} que nous avons explicitement spécifié dans sol[1].
Cette fonction d'interpolation n'a aucune expression analytique, mais nous pourrons quand même la nommer (disons vnum), la dessiner, la dériver, et l'approcher par une combinaison de fonctions de notre choix (grâce à la commande Fit).
Une fois avertis, nous passons à l'action :
sol[1]
vnum[s_]:= Evaluate[sol[1][[1,1,2]][s]]
graphevnum = Plot[Evaluate[vnum[s]],{s,-1,1}]
geod = ParametricPlot3D[x[s,vnum[s]]//Evaluate,{s,-1,1}]
geodsursphere = Show[sphere,geod,ViewPoint->{2,1.5,1.5},AxesLabel->{x,y,z}]
Ajoutons le point choisi à la figure :
pt = x[0,Pi/4]
dessinpoint = Show[ Graphics3D[{PointSize[.025],Point[pt]},ViewPoint->{2,1.5,1.5},AxesLabel->{x,y,z}]]
figurefinale = Show[sphere,geod,dessinpoint,ViewPoint->{2,1.5,1.5},AxesLabel->{x,y,z}]
Nous laissons au lecteur le soin de d'étudier la marge d'erreur commise par cette méthode. En particulier, vérifier que la courbe geod qu'on pourrait écrire simplement
c[s_]:= x[s,vnum[s]]
est bien sur la sphère, que c'est "presque" un cercle, que sa courbure géodésique est "presque nulle", que sa torsion est "presque nulle", etc.

Retour à ma page personnelle



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