# Definite Integral Fails due to Runtime Error

Here is the code that I have been writing; it's meant to predict the moment of inertia of a perfect cylinder tilted and offset from its axis of rotation, among other things. However, the definite_integral command fails with an error thrown by the gcd function. I suspect the issue is due to Sage's inner workings; however, I do not have the time or know-how to fully debug this issue. Any help would be appreciated.

from sage.symbolic.integration.integral import definite_integral
from sage.calculus.integration import numerical_integral
dSpringStretchedFully=10.625
dSpringUnstretchedLength=7.875
dArmatureInches=8
dSpringAnchor=4.75
thetaArmatureDefault=95
dSpringPitch=dSpringUnstretchedLength/sin(thetaArmatureDefault*pi/180)*sin(pi/180*(180-thetaArmatureDefault-(180/pi*arcsin(dSpringAnchor*sin(thetaArmatureDefault*pi/180)/dSpringUnstretchedLength))))
fSpringPoundsPerInch=65.69
kNewtonsPerPound=4.44822
kMetersPerInch=0.0254
fSpringNewtonsPerInch=fSpringPoundsPerInch*kNewtonsPerPound
dArmatureMeters=dArmatureInches*kMetersPerInch
dHalfRopeInches=dArmatureInches*cos(pi/180*(thetaArmatureDefault-90))
dHalfRopeMeters=dHalfRopeInches*kMetersPerInch
fSpringNewtonsPerMeter=fSpringNewtonsPerInch/kMetersPerInch
kSpringDirection=1
rRopeInches=0.5
rRopeMeters=rRopeInches*kMetersPerInch
var('x')
var('dHalfRope')
var('thetaArmatureToRope')
thetaPointToAxis(x,thetaArmatureToRope)=thetaArmatureToRope-(180/pi*arctan((x-rRopeMeters)/dArmatureMeters))
dPointFromFarAxis(x,dHalfRope)=((dHalfRope^2)+(x-rRopeMeters)^2)^(1/2)
dPointFromNearAxis(x,dHalfRope,thetaArmatureToRope)=((dArmatureMeters^2)+((dPointFromFarAxis(x,dHalfRope))^2)-2*dPointFromFarAxis(x,dHalfRope)*dArmatureMeters*cos(pi/180*thetaPointToAxis(x,thetaArmatureToRope)))^(1/2)
IMultiplier(x)=2*(((rRopeMeters^2)-(x-rRopeMeters)^2)^(1/2))
ILine(x,dHalfRope,thetaArmatureToRope)=((dPointFromNearAxis(x,dHalfRope, thetaArmatureToRope))^2)*IMultiplier(x)
ISlice(dHalfRope,thetaArmatureToRope)=definite_integral(ILine(x,dHalfRope,thetaArmatureToRope),x,0,2*rRopeMeters)
ISlice(3,4)

edit retag close merge delete

Sort by » oldest newest most voted

Given the fact that you have numerical parameters, i guess you want to obtain a numerical answer, not a symbolic one (which would not be very meaningful). So, instead of defining ISlice as a symbolic formula, you can make it a classical Python function, and use numerical_integral as follows:

sage: ISlice = lambda a,b : numerical_integral(ILine(x,a,b),0,2*rRopeMeters)
sage: ISlice(3,4)
(0.0039653374286522155, 3.5934168930166043e-07)


So, unless Sage is wrong somewhere, the numerical value of ISlice(3,4) is 0.0039653374286522155, with a precision of 3.6e-07.

Note: actually, you should do the same for the other functions, so that you do not have to use symbolic variable, nor symbolic functions:

sage: dSpringStretchedFully=10.625
....: dSpringUnstretchedLength=7.875
....: dArmatureInches=8
....: dSpringAnchor=4.75
....: thetaArmatureDefault=95
....: dSpringPitch=dSpringUnstretchedLength/sin(thetaArmatureDefault*pi/180)*sin(pi/180*(180-thetaArmatureDefault-(180/pi*arcsin(dSpringAnchor*sin(thetaArmatureDefault*pi/180)/dSpringUnstretchedLength))))
....: fSpringPoundsPerInch=65.69
....: kNewtonsPerPound=4.44822
....: kMetersPerInch=0.0254
....: fSpringNewtonsPerInch=fSpringPoundsPerInch*kNewtonsPerPound
....: dArmatureMeters=dArmatureInches*kMetersPerInch
....: dHalfRopeInches=dArmatureInches*cos(pi/180*(thetaArmatureDefault-90))
....: dHalfRopeMeters=dHalfRopeInches*kMetersPerInch
....: fSpringNewtonsPerMeter=fSpringNewtonsPerInch/kMetersPerInch
....: kSpringDirection=1
....: rRopeInches=0.5
....: rRopeMeters=rRopeInches*kMetersPerInch
....: thetaPointToAxis = lambda x,thetaArmatureToRope :thetaArmatureToRope-(180/pi*arctan((x-rRopeMeters)/dArmatureMeters))
....: dPointFromFarAxis = lambda x, dHalfRope : ((dHalfRope^2)+(x-rRopeMeters)^2)^(1/2)
....: dPointFromNearAxis = lambda x, dHalfRope,thetaArmatureToRope : ((dArmatureMeters^2)+((dPointFromFarAxis(x,dHalfRope))^2)-2*dPointFromFarAxis(x,dHalfRope)*dArmatureMeters*cos(pi/180*thetaPointToAxis(x,thetaArmatureToRope)))^(1/2)
....: IMultiplier = lambda x : 2*(((rRopeMeters^2)-(x-rRopeMeters)^2)^(1/2))
....: ILine = lambda x,dHalfRope,thetaArmatureToRope : ((dPointFromNearAxis(x,dHalfRope, thetaArmatureToRope))^2)*IMultiplier(x)
....: ISlice = lambda dHalfRope,thetaArmatureToRope : numerical_integral(ILine(x,dHalfRope,thetaArmatureToRope),0,2*rRopeMeters)
....: ISlice(3,4)
....:
(0.0039653374286522155, 3.5934168930166043e-07)

more

Thanks, that would definitely help. Three questions: the first output is the estimated value, and the second one is the uncertainty, correct? And secondly, do you have any idea why the symbolic integration fails? I'm asking because ISlice is still a few steps and integrations removed from the final answer that I seek, and keeping things symbolic until the final answer is calculated would probably result in better accuracy. Finally, if you don't really know why Sage is acting up, how would I integrate ISlice over a definite domain? Would Sage know how to take the resultant value automatically, or is there some extra step involved?

( 2017-11-01 23:34:54 +0200 )edit