Monte_carlo Integration and Gauss coil formula
Hi
When I try to compute the Gauss formula to calculate the number of turn's wraps between two solenoids.
if I don't put spaces between the turns (here a linear function, t*(turnSpace=0)
,
then the Integer value is inside the error range of the monte_carlo_integral()
function.
But if I introduce a space between the turns (turnSpace=2*pi/1024
),
then appears an error in the calculation of the Monte Carlo integral value which is out of the error range !.
(comment lines 61 and 62 to introduce a space between the turns of both coils)
Where is the error ? in my computation (most likely !)? in the monte_carlo_integral
function ?
with another software when integrating the function turnsNumberNumDerivative(t, s)
below between t = [0, 12*pi]
and s = [0, 14*pi]
do we find the integer value $42=6*7$ +- Error with the integer value in the error range?
turnsNumberNumDerivative(t,s)=-((sqrt(2)*pi*s + 1024*sqrt(3)*cos(s) - 512*sqrt(2)*sin(s) - 2457.6 *cos(t) - 1024)*(pi*(sqrt(3)*(sqrt(2)*pi - 512*sqrt(2)*cos(s)) + 1024*sin(s)) + 1228.8 *(sqrt(2)*pi + 512*sqrt(2)*cos(s))*cos(t)) - (pi*(sqrt(3)*(sqrt(2)*pi*s - 512*sqrt(2)*sin(s)) - 1024*cos(s) + 2457.6 *sin(t)) + 1228.8 *(sqrt(2)*pi*s - 2*pi*t + 512*sqrt(2)*sin(s))*cos(t))*(sqrt(2)*pi - 512*sqrt(2)*cos(s) - 1024*sqrt(3)*sin(s)) - 1228.8 *((sqrt(2)*pi + 512*sqrt(2)*cos(s))*(sqrt(3)*(sqrt(2)*pi*s - 512*sqrt(2)*sin(s)) - 1024*cos(s) + 2457.6 *sin(t)) - (sqrt(2)*pi*s - 2*pi*t + 512*sqrt(2)*sin(s))*(sqrt(3)*(sqrt(2)*pi - 512*sqrt(2)*cos(s)) + 1024*sin(s)))*sin(t))/(pi*(4*(sqrt(2)*pi*s - 2*pi*t + 512*sqrt(2)*sin(s))^2 + (sqrt(2)*pi*s + 1024*sqrt(3)*cos(s) - 512*sqrt(2)*sin(s) - 2457.6 *cos(t) - 1024)^2 + (sqrt(3)*(sqrt(2)*pi*s - 512*sqrt(2)*sin(s)) - 1024*cos(s) + 2457.6 *sin(t))^2)^(3/2))
new attempt to code button is now ok ;-), it seems I had to supress 00000000
print("https://www.maths.ed.ac.uk/~v1ranick/papers/ricca.pdf")
allVars=var('alpha,beta,gamma,'+
't,t0,t1,nt,shiftXt,shiftYt,shiftZt,Rt,turnSpaceT,'+
's,s0,s1,ns,shiftXs,shiftYs,shiftZs,Rs,turnSpaceS',domain='real')
assume(nt>0) ;assume(ns >0);
precis=50
#########################################
def addAxis(plt,textToPlt,xBound,yBound,zBound) :
# add axis for real part
plt += arrow3d((0, 0, 0), (xBound, 0, 0), color='blue')
plt += text3d("X", (xBound, 0, 0), color='blue')
plt += arrow3d((0, 0, 0), (0, yBound, 0), color='green')
plt += text3d("Y", (0, yBound, 0), color='green')
plt += arrow3d((0, 0, 0), (0, 0, zBound), color='red')
plt += text3d("Z", (0,0, zBound), color='red')
plt+=textToPlt
return plt
def RotateVector(V,anglesDic) :
#rotation angle alpha around X axis
RotalphaX=matrix(SR,[[1,0,0],[0,cos(alpha),-sin(alpha)],[0,sin(alpha),cos(alpha)]]).subs(anglesDic)
#rotation angle beta around Y axis
RotbetaY=matrix(SR,[[cos(beta),0,sin(beta)],[0,1,0],[-sin(beta),0,cos(beta)]]).subs(anglesDic)
#rotation angle gamma around Z axis
RotgammaZ=matrix(SR,[[cos(gamma),-sin(gamma),0],[sin(gamma),cos(gamma),0],[0,0,1]]).subs(anglesDic)
Vx= (RotalphaX*V)
Vxy= (RotbetaY*Vx)
Vxyz= (RotgammaZ*Vxy)
return vector(Vxyz)
anglesDict={alpha:0,beta:0,gamma:0}
anglesDics={alpha:pi/4,beta:0,gamma:pi/6}
#########################################
xyBound=1.2 ;
thck=8
gammat0 = vector(SR,[shiftXt+Rt*cos(t),shiftYt+Rt*sin(t),t*turnSpaceT]) # in xy plane
gammat=RotateVector(gammat0,anglesDict)
gammas0 = vector(SR,[shiftXs+Rs*cos(s),shiftYs+Rs*sin(s),s*turnSpaceS])
gammas= RotateVector(gammas0,anglesDics)
show("gammat : \t " ,gammat)
show("gammas : \t " ,gammas)
# set here below the number n of coil turns of the blue solenoid)
# the radius Rt and the turns space turnSpaceT
tDic={nt:6 , Rt:1 , turnSpaceT:2*pi/512 , shiftXt:1/2 , shiftYt:0 , shiftZt:0}
# change here below the number n of coil turns of the black solenoid)
# the radius Rs and the turns space turnSpaceS
sDic={ns:7 , Rs:1 , turnSpaceS:2*pi/512 , shiftXs:0 , shiftYs:0 ,shiftZs:0}
# no space between spires then Monte Carlo integration OK
#tDic={nt:6 , Rt:1 , turnSpaceT:0 , shiftXt:1/2 , shiftYt:0 , shiftZt:0}
#sDic={ns:7 , Rs:1 , turnSpaceS:0 , shiftXs:0 , shiftYs:0 ,shiftZs:0}
show("gammat with coil's values : \t " ,gammat.subs(tDic))
show("gammas with coil's values : \t " ,gammas.subs(sDic))
# set the circles centers
centerGammat=RotateVector(vector([shiftXt,shiftYt,shiftZt]).subs(tDic),anglesDict)
centerGammas=RotateVector(vector([shiftXs,shiftYs,shiftZs]).subs(sDic),anglesDics)
textToPlt = text3d("Gauss turns number", (xyBound, xyBound, xyBound), fontsize='300%', fontweight=800)
plt=parametric_plot3d(gammat.subs(tDic), (0,nt.subs(tDic)*2*pi),color="blue",thickness=thck)
plt+=parametric_plot3d(gammas.subs(sDic), (0,ns.subs(sDic)*2*pi),color="grey",thickness=thck)
plt+=addAxis(plt,textToPlt,1,1,1)
plt+=plot(centerGammat,color="blue")
plt+=point3d(centerGammat,size=20,color="blue")
plt+=plot(centerGammas,color="black")
plt+=point3d(centerGammas,size=20,color="black")
plt.show(viewer='threejs')
integrationFactorFunction=(1/(4*pi*((gammas-gammat).norm())^3)).subs(tDic).subs(sDic)
# here choose the integration function factor (integrationFactorFunction0 or integrationFactorFunction1)
M=matrix([ [ gammat.diff(t)[0],gammat.diff(t)[1],gammat.diff(t)[2] ],
[ gammas.diff(s)[0],gammas.diff(s)[1],gammas.diff(s)[2] ],
[(gammat[0]-gammas[0]), (gammat[1]-gammas[1]),(gammat[2]-gammas[2]) ] ])
precis=40 # precision approximation
# solenoid with space turn length
show(LatexExpr(r"tLen=\int^{b}_{a} \sqrt{ x_t′(t)^2 + y_t′(t)^2 + z_t′(t)^2 } dt \, \, \," +
r"sLen=\int^{b}_{a} \sqrt{ x_s′(t)^2 + y_s′(t)^2 + z_s′(t)^2 } dt "))
show(" ################ compute solenoid length ########################")
tLenDt=sqrt((M[0][0]^2 + M[0][1]^2 + M[0][2]^2).subs(tDic))
show(tLenDt)
solenoidTLen=integrate(tLenDt, (t, 0, nt.subs(tDic)*2*pi ) )
show("solenoid T length : \t ",solenoidTLen, " \t numerical value : \t ",solenoidTLen.n(precis))
sLenDt=sqrt((M[1][0]^2 + M[1][1]^2 + M[1][2]^2).subs(sDic))
show(sLenDt)
solenoidSLen=integrate(sLenDt, (s, 0, ns.subs(sDic)*2*pi ) )
show("solenoid S length : \t ",solenoidSLen, " \t numerical value : \t ",solenoidSLen.n(precis))
show("matrix M :")
show(M)
turnsNumber=integrationFactorFunction * M.det()
show("##################### turnsNumber #####################")
turnsNumberNumDerivative=(turnsNumber.subs(tDic)).subs(sDic)
textToPlts = text3d("surface to sum", (xyBound, xyBound, xyBound), fontsize='300%', fontweight=800,transparency=0.5)
var('x,y')
show("turns number fonction to integrate :",turnsNumberNumDerivative)
show(turnsNumberNumDerivative)
show("rough monte_carlo integral approximation : \t",
monte_carlo_integral(turnsNumberNumDerivative(x=t,y=s), [0,0], [ nt.subs(tDic)*2*pi, ns.subs(sDic)*2*pi], 10^8))
#for algo in ['plain', 'miser', 'vegas']: # abs tol 0.01
# show("rough monte_carlo integral approximation : \t",algo,
# monte_carlo_integral(turnsNumberNumDerivative(x=t,y=s), [0,0], [ nt.subs(tDic)*2*pi, ns.subs(sDic)*2*pi],
# 10^6, algorithm=algo))
show(r"""When I try to compute the Gaussian formula to calculate the number of turn's wraps between two solenoids.
if I don't put spaces between the turns (here a linear function,t*(turnSpace=0),
then the Integer value is inside the error range of the monte_carlo_integral() function.
But when I introduce a space between the turns(turnSpace=2*pi/1024), then appears an error in the computation
of the Monte Carlo integral value which is out of the error range !.
(comment lines 61 and 62 to introduce a space between the turns of both coils)
Where is the error ? in my computation ? in the monte_carlo_integral function ? in the Gauss formula ?""")
Please try to isolate the error in a minimal code, then post this minimal code. The question refers to the "intervals"
t
ands
, but when the code stars, these are variables. "Introducing a space" - and i suppose this means giving a non-zero value to the variableturnSpace
- is hard to be traced back in the code. Without an IDE it is even hard to digest the code, which is not PEP8-conformal, not documented for a neutral reader, and not structured so that the question would point to one point in it. Please help us with an isolated error, or just try to run the code in an IDE like eclipse, pycharm, spyder...Hi @dan_fulea ok it's hard to isolate the problem ! because maybe the trouble is in my calculation !
I added the function to integrate turnsNumberNumDerivative(t,s) explicitly at the end in the main body of tthe question, because it doesn't fit in comments.
with another software when integrating the function turnsNumberNumDerivative between t=[0,12pi] and s=[0,14pi]
do we find the integer value 42=6*7 +- Error with the integer value in the error range?
I don't understand when I select the function turnsNumberNumDerivative and click on code button the * disappear, the code is not accepted !. to get the function : in SageCell type: print(turnsNumberNumDerivative)