Ask Your Question
0

Monte_carlo Integration and Gauss coil formula

asked 2022-05-24 18:45:58 +0100

ortollj gravatar image

updated 2022-05-25 08:29:43 +0100

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 ?""")
edit retag flag offensive close merge delete

Comments

Please try to isolate the error in a minimal code, then post this minimal code. The question refers to the "intervals" t and s, but when the code stars, these are variables. "Introducing a space" - and i suppose this means giving a non-zero value to the variable turnSpace - 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...

dan_fulea gravatar imagedan_fulea ( 2022-05-24 19:09:44 +0100 )edit

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?

ortollj gravatar imageortollj ( 2022-05-24 19:40:38 +0100 )edit

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)

ortollj gravatar imageortollj ( 2022-05-24 19:52:47 +0100 )edit

1 Answer

Sort by » oldest newest most voted
0

answered 2022-05-26 08:42:27 +0100

ortollj gravatar image

I just understood, in fact there is no error! it is because when I put space between the turns, then my curve is not closed I do not return to the beginning of my curve, and therefore I do not respect the conditions of Gauss's theorem.

image description

edit flag offensive delete link more

Comments

Sorry for that question which was not a SageMath question but a simple math question ;-(

ortollj gravatar imageortollj ( 2022-05-26 08:56:17 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2022-05-24 18:45:58 +0100

Seen: 260 times

Last updated: May 26 '22