托马斯微积分 从入门到失望

import math
radius = 4.000
diameter = radius *2
# sphere r=4 R=8
# this spere is y=sqrt(16-x*x)
# per cylinder volume is PI*r*r*dtx
def sphere_function(xpos,dtx):
    return math.sqrt(16.00-xpos*xpos) * math.sqrt(16.00-xpos*xpos)*dtx

# @n is the is this sphere that will split to n piece along xpos
# if sphere r=4, slice is 8, dtx= R/8 = 1
def calculateVolume(n):
    v = 0
    dtx = diameter/n
    for j in range(0,n,1):
        xp = j*(diameter/n) - radius   #if sphere radius is 4,left plot is -4 ,right plot is 4
        v += sphere_function(xp,dtx)
    return v

# sphere volume use base function : 3/4 * PI * (R*R*R)
def ruleCacluateVolume(r):
    return 4/3.000 * r*r*r

if __name__ == "__main__":
    #split a sphere to 20 cylinder
    print "use 4   slice :"  ,calculateVolume(4)
    print "use 8   slice :"  ,calculateVolume(8)
    print "use 20  slice :"  ,calculateVolume(20)
    print "use 120 slice :"  ,calculateVolume(120)
    print "use 200 slice :"  ,calculateVolume(200)

import math

#
# curve function is y=sqrt(x)
# x range->0-5
#

maxRange = 5.0

def clinder_volume(xpos,dtx):
    return math.sqrt(xpos*xpos) * dtx

def curve_volume(n):
    v = 0.0
    dtx = maxRange/n
    for x in range(0,n,1):
        xpos = x*(maxRange/n)
        v += clinder_volume(xpos,dtx)
    return v

if __name__ == "__main__":
    print curve_volume(15)

# Trapezoidal
# S = 1/2(y0+ 2y1 + 2y2 + 2y3+...+ 2yn-1 + yn)

def Trapezoidal(down,up,n,func):
    if up==down:
        return 0.0
    h = float(up-down) / float(n)
    start = func(down)
    end = func(up)

    process = 0.0
    for dt in xrange(0,n+1,1):
        if dt == 0 or dt == n:
            continue
        process += 2 * func(down + dt * h)
    sum =  (start + end + process) * (h/2.0)
    return sum



# Simpson
# S = h/3(y0 + 4y1 + 2y2 + 4y3 + 2y4 + ... + 2yn-1 + yn)
# func is f(x)
def Simpson(down,up,n,func):
    if up==down:
        return 0.0
    h = float(up-down) / float(n)
    start = func(down)
    end = func(up)
    process = 0.0
    for dt in xrange(0,n+1,1):
        if dt == 0 or dt == n:
            continue
        # select the 1 3 5 7 9... index
        if dt%2 == 1:
            process += 4 * func(down + dt * h)
        # select the 2 4 6 8 10... index
        if dt%2 == 0:
            process += 2 * func(down + dt * h)

    sum =  (start + end + process) * (h/3.0)
    return sum



if __name__ == "__main__":

    # part1
    # fx = 5x^4 [0,2]  n=4
    func = lambda x:5*x*x*x*x
    T = Simpson(0,2,4,func)
    print T
    # part2
    # fx = x [1,2] n=4
    func2 = lambda x:x
    T2 = Simpson(1,2,4,func2)
    print T2

    # part3
    # fx = x*x
    func3 = lambda x:x*x
    T3 = Trapezoidal(1,2,4,func3)
    print T3

    func4 = lambda x:x*x + 1
    T_T4 = Trapezoidal(-1,1,4,func4)
    S_T4 = Simpson(-1,1,4,func4)
    print T_T4,S_T4

float ptsx[];
float ptsy[];
float ptsz[];
int npt = npoints(0);
resize(ptsx,npt);
resize(ptsy,npt);
resize(ptsz,npt);

for(int i=0;i
      
    
      
      
      
      
      
       
       
       
       
       
     
       )
{
    vector pos = point(
       
     
       
       
       
       
       0,
       
     
       
       
       
       
       "
       
     
       
       
       
       
       P
       
     
       
       
       
       
       "
       
     
       
       
       
       
       ,i); ptsx[i] =
       
     
       
       
       
       
        pos.x; ptsy[i] =
       
     
       
       
       
       
        pos.y; ptsz[i] =
       
     
       
       
       
       
        pos.z; } 
       
     
       
       
       
       
       float maxy =
       
     
       
       
       
       
        max(ptsy); 
       
     
       
       
       
       
       float miny =
       
     
       
       
       
       
        min(ptsy); printf(
       
     
       
       
       
       
       "
       
     
       
       
       
       
       %f,%f\n
       
     
       
       
       
       
       "
       
     
       
       
       
       
       ,miny,maxy); 
       
     
       
       
       
       
       float dttop = pow(maxy,
       
     
       
       
       
       
       3)  /  
       
     
       
       
       
       
       3.0
       
     
       
       
       
       
       ; 
       
     
       
       
       
       
       float dtbottom = pow(miny,
       
     
       
       
       
       
       3) / 
       
     
       
       
       
       
       3.0
       
     
       
       
       
       
       ; 
       
     
       
       
       
       
       float volume = dttop -
       
     
       
       
       
       
        dtbottom; adddetailattrib(geoself(),
       
     
       
       
       
       
       "
       
     
       
       
       
       
       cvolume
       
     
       
       
       
       
       ",volume);