Bonjour
j'ai besoin de calculer une integrale de surface un peu compliquée.
Aussi j'ai commencé par reviser les techniques d'integration.
et j'ai commencé par une integrale de volume (que je pensais simple !)
je me suis imposé un petit exercice a savoir calculer d'abord en coordonnées spheriques
une calotte sur une sphere et ensuite recalculer cette calotte quand elle est penchée de 45 degrés vers l'axe des y,
mais cette fois en coordonnées cartesiennes, mais ca ne colle pas et je ne trouve pas mon erreur.
question : est ce que le choix de mes borne d'integration est correcte ?
le code SageMath sur SageCell :
ci dessous le code sageMath
import matplotlib.pyplot as plt
import numpy as np
from sage.plot.line import Line
var('x,y,z,r,t,R,alpha,beta,gamma,theta,phi,h_up,h_low',domain='real')
assume(R>0)
assume (r>=0)
assume(r<=R)
Rnum=1
h_upNum=Rnum ;h_lowNum=sqrt(1/2)
numL=[R==Rnum,h_up==h_upNum,h_low==h_lowNum]
S = matrix(SR,[R*sin(theta)*cos(phi),R*sin(theta)*sin(phi),R*cos(theta)])
Snum=S.subs(numL)
SnumL=Snum.list()
xmin=-Rnum;xmax=Rnum;ymin=-Rnum;ymax=Rnum;zmin=-Rnum;zmax=Rnum
Plt=parametric_plot3d(SnumL, (theta,0,pi), (phi,0,2*pi), frame=False, color="blue",opacity=0.2)
f_h_up(x,y)=h_up
f_h_low(x,y)=h_low
f_h_up_num(x,y)=h_up.subs(numL)
f_h_low_num(x,y)=h_low.subs(numL)
# plot the three axis
Plt+=arrow((0,0,0), (Rnum,0,0), color='black')
Plt+=text3d("X",vector([Rnum,0,0])*1.1, color='black')
Plt+=arrow((0,0,0), (0,Rnum,0), color='brown')
Plt+=text3d("Y",vector([0,Rnum,0])*1.1, color='brown')
Plt+=arrow((0,0,0), (0,0,Rnum), color='red')
Plt+=text3d("Z",vector([0,0,Rnum])*1.1, color='red')
eqS=0==(R^2-x^2-y^2-z^2)
eqSt=eqS.subs(x=r*cos(t),y=r*sin(t)).simplify_trig().add_to_both_sides(r^2)
show(LatexExpr(r" R^2-x^2-y^2-z^2"),"\t passing in cylindrical coordinates : \t",eqSt)
#Cs(r,R)=(pi*r^2).subs(eqSt)
Cs=(pi*r^2).subs(eqSt)
show("circle surface :\t",LatexExpr(r" \,\pi \, R^2 "),"\t : Cs :\t",Cs)
Vs=integrate(Cs,z,h_low,h_up)
show("########################check for total sphere #####################################")
h_upNum=Rnum ;h_lowNum=-Rnum
numL=[R==Rnum,h_up==h_upNum,h_low==h_lowNum]
Plt+=plot3d(f_h_up_num,(x,xmin,xmax),(y,ymin,ymax), frame=False,color='green',opacity=0.2)
Plt+=plot3d(f_h_low_num,(x,xmin,xmax),(y,ymin,ymax), frame=False,color='pink',opacity=0.2)
zLow=-R
zUp=R
VsNum=Vs.subs(numL)
latexExprList=[]
latexExprList.append(LatexExpr(r"\text{vol of sphere} \,=\,"
+r"\int_{"+latex(zLow)+"}^{"+latex(zUp)+"}{"+ latex(Cs)+"} " \
+r"dz "))
show(latexExprList[-1])
show(" sphere total volume f(h_up,h_low) \t: \t",Vs)
show(" sphere total volume h_up=R and h_low==-R \t: \t",Vs.subs(h_low==-R,h_up=R),' \t numerical : \t',VsNum.n())
show("#######################check for sphere cap ###############################")
h_upNum=Rnum ;h_lowNum=1/sqrt(2)
numL=[R==Rnum,h_up==h_upNum,h_low==h_lowNum]
Plt+=plot3d(f_h_up_num,(x,xmin,xmax),(y,ymin,ymax), frame=False,color='green',opacity=0.2)
Plt+=plot3d(f_h_low_num,(x,xmin,xmax),(y,ymin,ymax), frame=False,color='pink',opacity=0.2)
zLow=h_low
zUp=R
VsNum=Vs.subs(numL)
latexExprList.append(LatexExpr(r"\text{vol of Cap only} \,=\,"
+r"\int_{"+latex(zLow)+"}^{"+latex(zUp)+"}{"+ latex(Cs)+"} " \
+r"dz "))
show(latexExprList[-1])
show("volume de la calotte Vs h_up=R and h_low==1/sqrt(2) \t : \t")
show(Vs,' \t numerical : \t',VsNum.n())
show(Plt,xmax=Rnum+0.1,ymax=Rnum+0.1)
S = matrix(SR,[R*sin(theta)*cos(phi),R*sin(theta)*sin(phi),R*cos(theta)])
Snum=S.subs(numL)
SnumL=Snum.list()
xmin=-Rnum;xmax=Rnum;ymin=-Rnum;ymax=Rnum;zmin=-Rnum;zmax=Rnum
Plt=parametric_plot3d(SnumL, (theta,0,pi), (phi,0,2*pi), frame=False, color="blue",opacity=0.2)
# cutting plan equation
p(x,y)=1-y
Plt+=plot3d(p,(x,xmin,xmax),(y,ymin,ymax), frame=False,color='green',opacity=0.3)
# plot the three axis
Plt+=arrow((0,0,0), (1,0,0), color='black')
Plt+=text3d("X",vector([1,0,0])*1.1, color='black')
Plt+=arrow((0,0,0), (0,1,0), color='brown')
Plt+=text3d("Y",vector([0,1,0])*1.1, color='brown')
Plt+=arrow((0,0,0), (0,0,1), color='red')
Plt+=text3d("Z",vector([0,0,1])*1.1, color='red')
assume((1-y)>=0)
assume((R^2-x^2-y^2)>=0)
assume(2*R^2-1>0)
eqS=0==(R^2-x^2-y^2-z^2)
eqSt=eqS.subs(x=r*cos(t),y=r*sin(t)).simplify_trig().add_to_both_sides(r^2)
show(LatexExpr(r" R^2-x^2-y^2-z^2"),"\t pass in cylindrical coordinates : \t",eqSt)
eqStt=solve(eqS,z^2)[0].rhs()
eqSttt=eqStt==(1-y)^2
eqSx=eqSttt.add_to_both_sides(-R^2+y^2).multiply_both_sides(-1)
Sx=solve(eqSx,x)
show("eqSttt : ",eqSttt,"\t Ellipse : ",LatexExpr(r"\frac{(x - u)^2}{a^2} + \frac{(y - v)^2}{b^2} = 1"))
show("eqSx : ",eqSx)
# we choose to integrate under the plane 1-y, on upward semi-sphere and y >= 0
#show(solve([,z^2==(1-y)^2],X,Y))
#show(integrate(x,integrate(integrate(r,r,0,R-y),y,0,1),-1,1))
#Plt=parametric_plot3d(SnumL, (theta,0,pi), (phi,0,2*pi), frame=False, color="blue",opacity=0.1)
#R=rNum
#eqE=1==x^2/R^2+(y-1/2)^2/(R^2/4)
eqE=solve(eqSttt,y)
show("Ellipse Equation : ",eqE)
eq0=eqE[1].multiply_both_sides(2).add_to_both_sides(-1)
show("eq0 : ",eq0)
spanX=solve(eq0.rhs()==0,x)
show("Ellipse x span :",spanX,"\t ok")
Plt+=plot(eqE[0].rhs().subs(numL),(x,spanX[0].rhs().subs(numL),spanX[1].rhs().subs(numL)),color='black')
Plt+=plot(eqE[1].rhs().subs(numL),(x,spanX[0].rhs().subs(numL),spanX[1].rhs().subs(numL)),color='black')
#Plt+=plot(spanX[0].rhs().subs(numL).n(),(x,-R,R),color='grey')
pts00=[[spanX[0].rhs().subs(numL),0] ,[spanX[0].rhs().subs(numL),Rnum]]
Plt+=line(pts00,color='blue')
pts01=[[spanX[1].rhs().subs(numL),0] ,[spanX[1].rhs().subs(numL),Rnum]]
Plt+=line(pts01,color='blue')
# volume above the ellipse in x ,y plan
#Plt+=list_plot([0,1/2,1/sqrt(2)],color='black')
vH=vector([0,1/2,1/2])
Plt+=plot(vH, color='pink')
Plt+=text3d("h",vH*1.5, color='pink',fontweight='bold', fontstyle='italic', fontsize=20,opacity=1)
from sage.plot.plot3d.shapes2 import Line
Plt+=Line([(0, 1, 0),(0, 0, 1)],color="pink" )
zUp=R
zLow=1-y
zZero=0
xLow= Sx[0].rhs()
show("x(y) for lower part of the ellipse : \t",Sx[0].rhs())
xUp= Sx[1].rhs()
show("x(y) for upper part of the ellipse : \t",Sx[1].rhs())
yLow=eqE[0].rhs()
yUp=eqE[1].rhs()
show('xLow : \t',xLow,'\t xUp : \t',xUp,'\t yLow : \t',yLow,'\t yUp : \t',yUp,'\t zLow : \t',zLow,'\t zUp : \t',zUp)
sphereVolume=integrate(integrate(integrate(r,r,0,sqrt(R^2-z^2)),z,-R,R),t,0,2*pi)
show("radius R sphereVolume : ",sphereVolume," \t num : ",sphereVolume.subs(numL).n())
show("########################################################################" )
zLow=0
zUp=R^2-x^2-y^2
yLow=0
yUp=1
latexExprList.append(LatexExpr(r"\text{vol total above ellipse in x,y plan} \,=\,"
+r"\int_{"+latex(zLow)+"}^{"+latex(zUp)+"} " \
+r"\int_{"+latex(xLow)+"}^{"+latex(xUp)+"} " \
+r"\int_{"+latex(yLow)+"}^{"+latex(yUp)+"} " \
+r"dz\, dx\, dy"))
show(latexExprList[-1])
volumeAboveEllipsoidalCylinder=integrate(integrate(integrate(1,z,zLow,zUp),x,xLow,xUp),y,yLow,yUp)
show("Volume total Above Ellipsoidal Cylinder in xy plan : ",volumeAboveEllipsoidalCylinder,"\t numerical : ",volumeAboveEllipsoidalCylinder.subs(numL).n())
show("########################################################################" )
zLow=1-y
latexExprList.append(LatexExpr(r"\text{vol of sphere Cap only} \,=\,"
+r"\int_{"+latex(zLow)+"}^{"+latex(zUp)+"} " \
+r"\int_{"+latex(xLow)+"}^{"+latex(xUp)+"} " \
+r"\int_{"+latex(yLow)+"}^{"+latex(yUp)+"} " \
+r"dz\, dx\, dy"))
show(latexExprList[-1])
volumeCapAbovePlan=integrate(integrate(integrate(1,z,zLow,zUp),x,xLow,xUp),y,yLow,yUp)
show("Volume Cap Above plan 1-y : ",volumeCapAbovePlan,"\t numerical : ",volumeCapAbovePlan.subs(numL).n())
show("########################################################################" )
zLow=0
zUp=1-y
latexExprList.append(LatexExpr(r"\text{vol between plan x,y and plan 1-y delimited by ellipse } \,=\,"
+r"\int_{"+latex(zLow)+"}^{"+latex(zUp)+"} " \
+r"\int_{"+latex(xLow)+"}^{"+latex(xUp)+"} " \
+r"\int_{"+latex(yLow)+"}^{"+latex(yUp)+"} " \
+r"dz\, dx\, dy"))
show(latexExprList[-1])
volumeBelowTruncatedCylinder=integrate(integrate(integrate(1,z,zLow,zUp),x,xLow,xUp),y,yLow,yUp)
show("Volume Truncated ellipsoid Cylinder between 1-y plan and x,y plan : ",volumeBelowTruncatedCylinder,"\t numerical : ",volumeBelowTruncatedCylinder.subs(numL).n())
show("########################################################################" )
show(Plt,xmax=Rnum+0.1,ymax=Rnum+0.1)
Visiblement tes bornes d'intégration ne vont pas : y varie de 0 à R et pour x ce n'est pas homogène.
mais y varie bien de zero a R, y doit bien parcourir tout le demi petit axe de l'ellipse ?
(pour les coordonnées cartesiennes, je suis sur de mon integrale en coordonnées polaires)
Mais regarde donc ce que tu as écrit ! Tu as écrit une intégrale de 0 à 1.
Et pour x, je le répète, tes bornes ne sont pas homogènes.
oops ok pour l'erreur du 1 (comme R a pour valeur 1 dans un premier temps , je me suis trompé et j'ai mis 1 a la place de R).
par contre je ne comprends pour l'instant pas ta remarque sur l'homogeneité des bornes pour x, pourrais tu preciser stp ?
je ne sais pas si il est possible de changer yUp=1 en yUp= R dans le code joint dans le premier post ?
x, y , R sont des longueurs. tes bornes d'intégration doivent être homogènes à des longueurs (au sens des dimensions de grandeurs en physique). Elles ne le sont pas du tout.
pourrais tu me confirmer que mes bornes pour y, et z sont ok maintenant ?
et qu'il me reste les bornes pour x a trouver ?
je ne suis pas d'accord avec toi mes bornes pour x sont correctes et correspondent bien a une dimension
avec le choix de mes bornes initiales pour x, je me retrouve avec une fonction:
qui n'est pas integrable algebriquement en y
mais qui l'est en numerique.
la preuve ci dessous:
f(y)=-(R^2 - y^2)*arcsin(-sqrt(R^2 - 2*y^2 + 2*y - 1)/sqrt(R^2 - y^2)) + sqrt(R^2 - 2*y^2 + 2*y - 1)*(2*y + sqrt(y^2 - 2*y + 1) - 2)
g(x)=f.subs(y==x).subs(R==1)
show(g)
numerical_integral(g, 0, 1, max_points=100)
merci pour ton aide GBMZ.
en fait tout a l'origine (pas celle que j'ai ecrit dans mon premier post)ma premiere ecriture etait bonne:
mais comme je n'arrivais pas a integrer,
j'avais fais un essai en supprimant la racine carré de la borne pour z,
et je ne me souvenait plus l'avoir enlever !!
Non, ça ne va toujours pas : tes bornes pour x et ta borne inférieure pour z ne sont pas homogènes à des longueurs. J'ai l'impression que dans tes calculs tu prends des fois R et des fois 1 : c'est typique pour le z=1-y, alors que ton dessin montre qu'il s'agit de z=R-y.
oui je suis d'accord pour le R-y
c'est du fait que je veux faire un calcul generique pour R mais dans ma tete je reviens toujour a R=1 !
voila:
Là au moins, il n'y a pas de problème d'homogénéité. Ça ne garantit pas la correction, mais ça rassure.
Coucou c'est re-moi !
(cest un peu la suite de: techniques integration volume
Je me pose des question sur une de mes questions
sur SageMath forum
(ca fait un peu boite de cammembert,d'ailleurs ca va continuer question sur ma question sur ma question !)
Personne n'a répondu ,je sais bien qu'il n'y a aucune obligation de répondre a une question mathematique,
et je n'exige pas de réponse, d'ailleurs peut etre que ma question mathematique n'est pas pertinente .
et donc je me pose des questions sur ma question,j'aurais voulu un avis
(on a le droit de cocher plusieur cases)
1er) je me suis planté dans mon calcul d'integrale ?
2em) ma question est loufoque ? du genre:
(d'ailleurs je soupconne un troll derriere cette question)
3em) ma question n'est pas interessante ?
4em) il n'y a pas de solution simple a ma question
5em)autre
j'explique ce que je voulais faire dans ce code (sans doute naivement et voué a l'echec,(Liouville)
mais ca c'est pas grave et ca n'est pas le sujet, c'est juste par curiosité, meme de voir ou ca coince ) :
je voulais tenter de d'obtenir une longueur algebrique d'un arc d'ellipse a l'aide du theoreme d'Harriot.
etape 1) verifier le code afin qu'on retrouve bien numeriquement l'aire de surface d'un triangle spherique quelconque
situé dans le 1er octant ,calculée par le theorem d'Harriot et par l'integral de surface dy dx
1er integrale dy algebrique, 2em dx numerique.
etape 2) me servir de la valeur algebrique obtenu par Harriot
pour tenter de retrouver une valeur algebrique d'un arc d'ellipse
je peux peut etre mettre ici le code SageMath complet ? mais il est un peu long.
et il y a un défaut on ne peut le resseter sans que ca plante le notebook.
*** message déplacé ***
j'ai honte je viens de m'apercevoir qu'on avait en fait répondu a ma question !!
dans la suite de mes commentaires il sufisait de cliquer sur more comments !!!
est il possible d'effacer ce post de ile aux maths?
*** message déplacé ***
Vous devez être membre accéder à ce service...
Pas encore inscrit ?
1 compte par personne, multi-compte interdit !
Ou identifiez-vous :