2447.sagews

# 薛鹏翔 20100915
# 定义给定高度height和半径radius的圆缺的面积
def HemiCircle(radius,height):
    var('r y h')
    assume(r>0)
    assume(0<h<=2*r)
    f=2*( r^2 - (y-r)^2 )^(1/2)
    area=f.integrate(y,0,h)
    if height>2*radius:
        height=2*radius
    if height<0:
        height=0
    return area(r=radius,h=height)

# 定义给定高度height和长半径longradius、短半径shortradius的椭圆面积
def HemiEllipse(longradius,shortradius,height):
    var('a b y h')
    assume(a>0)
    assume(b>0)
    assume(0<h<2*b)
    f=2*( ( 1 - (y-b)^2/(b^2) )*(a^2) )^(1/2)
    area=f.integrate(y,0,h)
    if height>2*shortradius:
        height=2*shortradius
    if height<0:
        height=0
    return area(a=longradius,b=shortradius,h=height)

# 定义给定油罐侧切面椭圆长半径longradius、短半径shortradius、显示油高height、油浮左侧长度leftlength、油浮右侧长度rightlength、和倾斜角度alpha时的体积
def TankVolume(longradius,shortradius,leftlength,rightlength,height,alpha):
    var('x h')
    assume(x>0)
    assume(h>0)
    f=HemiEllipse(longradius,shortradius,height-x*tan(alpha))
    if height<=(leftlength+rightlength)*tan(alpha):
        area=f.integrate(x,0,height/tan(alpha))
    else:
        area=f.integrate(x,0,leftlength+rightlength)
    return area(h=height+leftlength*tan(alpha))

# 定义球缺半径为radius、宽度为length、高度为height和倾斜角度为alpha的体积
def hemispherical(radius,length,height,alpha):
    var('x r h H')
    assume(0<x<radius-length)
    Height=2*(radius^2-(radius-length)^2)^(1/2)
    if height>Height:
        height=Height
    if height<0:
        height=0
    # 计算分割点
    if height==Height:
        xx=0
    if height==0:
        xx=0
    if 0<height<Height:
        var('r l h H a x')
        xx=(-1/2*(H*tan(a) - 2*h*tan(a) - 2*l + 2*r - sqrt(-4*l^2*tan(a)^2 +
8*h*l*tan(a) + 8*(l*tan(a)^2 - h*tan(a))*r - 4*(l*tan(a) - r*tan(a) -
h)*H - H^2 - 4*h^2 + 4*r^2))/(tan(a)^2 + 1))(r=radius,l=length,h=height,H=Height,a=alpha)#事先计算好该表达式
    #计算分割点左端积分
    assume(r>0)
    assume(0<h<2*r)
    f=HemiCircle(r,h)
    rr=(radius^2-(radius-length+x)^2)^(1/2)
    hh=height-Height/2-x*tan(alpha)+rr
    g=f(h=hh,r=rr)
    volume=integral_numerical(g,0,xx)[0]
    #volume=integrate(g,x,0,xx)
    #计算分割点右端积分
    g=f(h=2*rr,r=rr)
    volume+=integral_numerical(g,xx,length)[0]
    #volume+=integrate(g,x,xx,length)
    return volume

# 定义旋转alpha角度和beta角度后油罐的容量
def Tank2Volume(radius,sphericalradius,length,leftlength,rightlength,height,alpha):
    #中间部分
    volume=TankVolume(radius,radius,leftlength,rightlength,height,alpha)
    #左侧油罐的体积
    leftheight=height+leftlength*tan(alpha)
    volume+=hemispherical(sphericalradius,length,leftheight,-alpha)
    #右侧油罐的体积
    rightheight=height-rightlength*tan(alpha)
    volume+=hemispherical(sphericalradius,length,rightheight,alpha)
    return volume(h=height+leftlength*tan(alpha))

# 定义第一问的体积V和高度h、旋转角度alpha的函数关系
var('h alpha')
assume(tan(alpha/180*pi)>0)
assume(h<0.6)#V和V0的表达式中惠出现arcsin，SAGE要求明确，
V=TankVolume(1.78/2,1.2/2,0.4,2.05,h,alpha/180*pi)
V0=TankVolume(1.78/2,1.2/2,0.4,2.05,h,0)#零度时单独定义，原因为积分式子中包含tan的倒数

# 水平入油
hh=[0.15902,0.17614,0.19259,0.20850,0.22393,0.23897,0.25366,0.26804,0.28216,0.29603,0.30969,0.32315,0.33644,0.34957,0.36256,0.37542,0.38816,0.40079,0.41332,0.42576,0.43812,0.45040,0.46262,0.47478,0.48689,0.49895,0.51097,0.52295,0.53490,0.54682,0.55872,0.57061,0.58248,0.59435,0.60622,0.61809,0.62996,0.64185,0.65375,0.66567,0.67763,0.67854,0.69053,0.69082,0.70285,0.71491,0.72703,0.73919,0.75142,0.76370,0.76416,0.77653,0.78899,0.80154,0.81419,0.82695,0.83983,0.85284,0.86600,0.87932,0.89282,0.89284,0.90653,0.92045,0.93461,0.94905,0.96380,0.97891,0.99443,1.01043,1.02699,1.04425,1.06237,1.08159,1.10233,1.12532,1.15236,1.19349]
vv=[0.312,0.362,0.412,0.462,0.512,0.562,0.612,0.662,0.712,0.762,0.812,0.862,0.912,0.962,1.012,1.062,1.112,1.162,1.212,1.262,1.312,1.362,1.412,1.462,1.512,1.562,1.612,1.662,1.712,1.762,1.812,1.862,1.912,1.962,2.012,2.062,2.112,2.162,2.212,2.262,2.312,2.31583,2.36583,2.36706,2.41706,2.46706,2.51706,2.56706,2.61706,2.66698,2.66883,2.71883,2.76883,2.81883,2.86883,2.91883,2.96883,3.01883,3.06883,3.11883,3.16883,3.16891,3.21891,3.26891,3.31891,3.36891,3.41891,3.46891,3.51891,3.56891,3.61891,3.66891,3.71891,3.76891,3.81891,3.86891,3.91891,3.96891]
p=list_plot(zip(hh,vv))
p+=plot(V0(h=x),(x,0,1.2),color='red')
p.show(frame=true)

#绝对误差
verror=[(V0(h=hh[each])-vv[each]) for each in range(len(hh))]
p=list_plot(zip(hh,verror))
p.show(frame=true)

# alpha=0 beta=0时刻度对应的油量
print Tank2Volume(1.5,1.625,1,2,6,2632.23/1000,0).n()
print Tank2Volume(1.5,1.625,1,2,6,2624.3/1000,0).n()
print Tank2Volume(1.5,1.625,1,2,6,2620.67/1000,0).n()
print Tank2Volume(1.5,1.625,1,2,6,2610.29/1000,0).n()
print Tank2Volume(1.5,1.625,1,2,6,2606.61/1000,0).n()

60.4488956583381 60.3114390818059 60.2480411121023 60.0651234056804 59.9996994123829

#穷举法计算alpha和beta，对数据分析发现最大误差(0.6)出现在显示油高为1653.73附近，由此，控制该部分误差达到最小，则全局误差近似达到最小
for k in range(0,5):
    for m in range(0,3):
        alpha=k/180*pi
        beta=m/180*pi
        rotateheight=(2.63223-1.5)*cos(beta)+1.5
        tankorigin=Tank2Volume(1.5,1.625,1,2,6,rotateheight,n(alpha)).n()
        rotateheight=(1.65373-1.5)*cos(beta)+1.5
        tank=Tank2Volume(1.5,1.625,1,2,6,rotateheight,n(alpha)).n()
        print "角度alpha：",alpha/pi*180,"beta:",beta/pi*180,"误差",tankorigin-tank-24.42887

角度alpha： 0 beta: 0 误差 -0.658941779641172 角度alpha： 0 beta: 1 误差 -0.661256402618381 角度alpha： 0 beta: 2 误差 -0.668203653080003 角度alpha： 1 beta: 0 误差 -0.0595563704500854 角度alpha： 1 beta: 1 误差 -0.0620984346606441 角度alpha： 1 beta: 2 误差 -0.0697275159114064 角度alpha： 2 beta: 0 误差 0.416291608802741 角度alpha： 2 beta: 1 误差 0.413552961676249 角度alpha： 2 beta: 2 误差 0.405334542584772 角度alpha： 3 beta: 0 误差 0.775244380360387 角度alpha： 3 beta: 1 误差 0.772334044912956 角度alpha： 3 beta: 2 误差 0.763600921494774 角度alpha： 4 beta: 0 误差 1.01887859189684 角度alpha： 4 beta: 1 误差 1.01581752294704 角度alpha： 4 beta: 2 误差 1.00663252654408

#由上一次循环发现误差最小(0.059)出现在alpha=1,beta=0时，为提高精度，现步长选择为0.1
for k in range(5,15):
    for m in range(0,5):
        alpha=k/1800*pi
        beta=m/1800*pi
        rotateheight=(2.63223-1.5)*cos(beta)+1.5
        tankorigin=Tank2Volume(1.5,1.625,1,2,6,rotateheight,n(alpha)).n()
        rotateheight=(1.65373-1.5)*cos(beta)+1.5
        tank=Tank2Volume(1.5,1.625,1,2,6,rotateheight,n(alpha)).n()
        print "角度alpha：",alpha/pi*180,"beta:",beta/pi*180,"误差",tankorigin-tank-24.42887

角度alpha： 1/2 beta: 0 误差 -0.343187788841515 角度alpha： 1/2 beta: 1/10 误差 -0.343212113206551 角度alpha： 1/2 beta: 1/5 误差 -0.343285086613907 角度alpha： 1/2 beta: 3/10 误差 -0.343406710001528 角度alpha： 1/2 beta: 2/5 误差 -0.343576984930888 角度alpha： 3/5 beta: 0 误差 -0.283938691146961 角度alpha： 3/5 beta: 1/10 误差 -0.283963240869685 角度alpha： 3/5 beta: 1/5 误差 -0.284036890345323 角度alpha： 3/5 beta: 3/10 误差 -0.284159640496814 角度alpha： 3/5 beta: 2/5 误差 -0.284331492862041 角度alpha： 7/10 beta: 0 误差 -0.225961708526118 角度alpha： 7/10 beta: 1/10 误差 -0.225986480232120 角度alpha： 7/10 beta: 1/5 误差 -0.226060795652021 角度alpha： 7/10 beta: 3/10 误差 -0.226184655695011 角度alpha： 7/10 beta: 2/5 误差 -0.226358061874645 角度alpha： 4/5 beta: 0 误差 -0.169245867298446 角度alpha： 4/5 beta: 1/10 误差 -0.169270857700280 角度alpha： 4/5 beta: 1/5 误差 -0.169345829203802 角度alpha： 4/5 beta: 3/10 误差 -0.169470782703694 角度alpha： 4/5 beta: 2/5 误差 -0.169645719690788 角度alpha： 9/10 beta: 0 误差 -0.113780731783802 角度alpha： 9/10 beta: 1/10 误差 -0.113805937677011 角度alpha： 9/10 beta: 1/5 误差 -0.113881555650277 角度alpha： 9/10 beta: 3/10 误差 -0.114007586584389 角度alpha： 9/10 beta: 2/5 误差 -0.114184031947222 角度alpha： 1 beta: 0 误差 -0.0595563704500854 角度alpha： 1 beta: 1/10 误差 -0.0595817887073657 角度alpha： 1 beta: 1/5 误差 -0.0596580437683407 角度alpha： 1 beta: 3/10 误差 -0.0597851365005120 角度alpha： 1 beta: 2/5 误差 -0.0599630683492087 角度alpha： 11/10 beta: 0 误差 -0.00656332437099749 角度alpha： 11/10 beta: 1/10 误差 -0.00658895193779685 角度alpha： 11/10 beta: 1/5 误差 -0.00666583492315098 角度alpha： 11/10 beta: 3/10 误差 -0.00679397418112515 角度alpha： 11/10 beta: 2/5 误差 -0.00697337113507146 角度alpha： 6/5 beta: 0 误差 0.0452074241622462 角度alpha： 6/5 beta: 1/10 误差 0.0451815902717208 角度alpha： 6/5 beta: 1/5 误差 0.0451040883199347 角度alpha： 6/5 beta: 3/10 误差 0.0449749174656766 角度alpha： 6/5 beta: 2/5 误差 0.0447940763071735 角度alpha： 13/10 beta: 0 误差 0.0957644744734303 角度alpha： 13/10 beta: 1/10 误差 0.0957384371801631 角度alpha： 13/10 beta: 1/5 误差 0.0956603250244719 角度alpha： 13/10 beta: 3/10 误差 0.0955301371778923 角度alpha： 13/10 beta: 2/5 误差 0.0953478722601737 角度alpha： 7/5 beta: 0 误差 0.145115971898523 角度alpha： 7/5 beta: 1/10 误差 0.145089734059955 角度alpha： 7/5 beta: 1/5 误差 0.145011020272303 角度alpha： 7/5 beta: 3/10 误差 0.144879829719745 角度alpha： 7/5 beta: 2/5 误差 0.144696161042990

#对照上面的结果，发现答案为alpha=1.1,beta=0,误差为0.00656（此精度应达到要求了，当然可以再自1.1左右搜索更精确的解，这里就不累赘了）
#穷举法是一种最准确却最笨的方法，这里由于考虑到体积对参数alpha、beta的连续性，采用这样的多尺度分析（由粗至细）可以加快速度，不过可惜的是如上的两次计算仍旧消耗了近10分钟时间
#通过这次的题目，看到了Matlab的无奈，毕竟Matlab不是数学软件。呵呵，折腾并快乐着，从mathematica到matlab，再到Sage之路，不知未来如何（GSL或许是下一个选择）。
#通过最终对数据的计算，发现题目提供的数据为计算机模拟，即按照事先给定的公式角度算处的数据，非真实数据，第一问的相对误差恒定可能与产生数据时的计算精度和油罐尺寸精度油罐
#经过重新计算，发现题目当中数据应该是给错了，第一问长半径应该是1.72

#前边使用的计算分割点表达式
var('r l h H a x')
xx=solve(r^2-(r-l+x)^2==(h-H/2-x*tan(a))^2,x,solution_dict=True)
xx

[{x: -1/2*(H*tan(a) - 2*h*tan(a) - 2*l + 2*r + sqrt(-4*l^2*tan(a)^2 + 8*h*l*tan(a) + 8*(l*tan(a)^2 - h*tan(a))*r - 4*(l*tan(a) - r*tan(a) - h)*H - H^2 - 4*h^2 + 4*r^2))/(tan(a)^2 + 1)}, {x: -1/2*(H*tan(a) - 2*h*tan(a) - 2*l + 2*r - sqrt(-4*l^2*tan(a)^2 + 8*h*l*tan(a) + 8*(l*tan(a)^2 - h*tan(a))*r - 4*(l*tan(a) - r*tan(a) - h)*H - H^2 - 4*h^2 + 4*r^2))/(tan(a)^2 + 1)}]