Monte_carlo Integration and Gauss coil formula

edit

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?

code on SageCell

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 ?""")