Definite Integral Fails due to Runtime Error

asked 2017-11-01 01:11:46 +0100

updated 2017-11-01 19:02:08 +0100

tmonteil
27423 ●31 ●204 ●517 http://wiki.sagemath.o...

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
fSpringMinimumLoad=0
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)

answered 2017-11-01 21:42:11 +0100

tmonteil
27423 ●31 ●204 ●517 http://wiki.sagemath.o...

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
....: fSpringMinimumLoad=0
....: 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)

edit flag offensive delete link

Comments

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?

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

add a comment

Your Answer

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

Definite Integral Fails due to Runtime Error

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Definite Integral Fails due to Runtime Error edit

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Definite Integral Fails due to Runtime Error