Ask Your Question
1

Definite Integral Fails due to Runtime Error

asked 2017-10-31 19:11:46 -0600

Benji_Zhang gravatar image

updated 2017-11-01 13:02:08 -0600

tmonteil gravatar image

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

1 answer

Sort by ยป oldest newest most voted
1

answered 2017-11-01 15:42:11 -0600

tmonteil gravatar image

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 more

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 gravatar imageBenji_Zhang ( 2017-11-01 17:34:54 -0600 )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: 2017-10-31 11:17:36 -0600

Seen: 33 times

Last updated: Nov 01