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é ***
merci de l'avoir déplacé et mis a la suite de integral de volume !
deja parceque c'est effectivement la suite et de plus c'est plus discret
Bonjour
le sujet maintenant est integral curviligne.
je croyais avoir compris , ce qui etait calculé dans la video :,, j'ai fais un commentaire sur la video (voir le commentaire de Zozio Descimes
ortollj et Zozio Descimes ne font qu'un!) mais maintenant je suis perdu. que represente la quantité obtenu dans le calcul integrale quand une surface rebique comme ci dessous et que le segment du chemin d'integration traverse la surface !!.
autrement dit que representent les deux quantités 13.766 (bas du segment jusqu'au point rouge), -13.766 (point rouge au haut du segment) ?
code on SageCell:
Ta représentation graphique ne va pas du tout : tu représentes la surface d'équation au lieu de représenter la surface .
Tu n'es pas d'accord avec quoi ? La surface représentée en bleu sur ta figure 3d n'est pas la surface d'équation ???
Enfin, comme on a pu le voir dans la discussion précédente, tu commences par dire que tu n'es pas d'accord, pour finalement te rendre compte que j'ai raison. Il suffit que j'attende un petit peu.
En fait je me suis fais un noeud au cerveau pour pas grand chose.
Je pense apres reflexion qu'il ne s'agit pas de hasard, mais que les 2 bornes [[3, 6, 0], [-1, 1, 0]],
avaient du originalement etre choisies pour que l'integrale soit egale a zero.
J'ai rajouté au code apres la ligne 100 les lignes suivantes:
lineIntegralLowTogreenPt=integral(lineIntegrand,t,numL[0].rhs(),St[0][0].rhs())
show("lineIntegral low point to green Pt : ",lineIntegralLowTogreenPt," \t numerical : ",
round(lineIntegralLowTogreenPt,7))
lineIntegralGreenToRedPt=integral(lineIntegrand,t,St[0][0].rhs(),St[1][0].rhs())
show("lineIntegral green point to red Pt : ",lineIntegralGreenToRedPt," \t numerical : ",
round(lineIntegralGreenToRedPt,7))
Bon il y a encore un truc qui me chiffone, c'est le signe des surfaces !!
ils devraient etre inversés.
J'ai fait en SageMath, puisque c'est ce que tu utilises. C'est tout bête :
x,y,t=var('x,y,t')
surf=plot3d(3*x^2-2*y,(x,-1,3),(y,1,6),opacity=.7)
seg=line([(-1,1,0),(3,6,0)],thickness=5, color='green')
cour=parametric_plot3d([-1+4*t,1+5*t, 3*(4*t-1)^2-2*(5*t+1)],\
(t,0,1),thickness=5, color='red')
show(surf+seg+cour, aspect_ratio=[3,3,1])
Vous devez être membre accéder à ce service...
Pas encore inscrit ?
1 compte par personne, multi-compte interdit !
Ou identifiez-vous :