nodes - 复合梯形规则和标准正态分布所需的节点数

Question

我正在使用不同的数值方法来近似正态分布。我的第一种方法是梯形规则。为什么我的节点数量需要表现得如此奇怪？我使用了一个表格来找到 n 的最小值，这样误差就会小于 0.0001。我会期待一个模式。有没有，我错过了？标准差为 1 的 n 为 6。标准差 2 的 n 为 54。标准差 3 的 n 为 29。

我的代码如下：

s=1---------------------------------------------------------
print("For Standard Deviation, s, of 1, we have a=-1 and b=1")


def CompositeTrapezoidal(f, a, b, n):
    #compute the size of each subinterval
    h = (b-a)/n
    #here we add all the terms
    sum = 1/2*f(a)
    #add the values for all intermediate nodes
    for i in [1,2,..,n-1] :
        sum += f(a+i*h)
    #add the last node
    sum +=1/2*f(b)
    #multiply by h to complete the computation
    sum = sum*h
    #return the sum
    return sum
#define the function we will be using
f(x) =(e^-((x^2)/2))/(sqrt(2*pi))
#print the top part of the table
print ("n", "\t","trapezoidal" , "\t\t", "error")
print ("--------------------------------------------------")
#call the function several times
for power in [11,12,..,45]:
    n = power
    trapezoidal = CompositeTrapezoidal(f, -1, 1, n).n()
    error = (trapezoidal - integral(f(x), x,-1,1)).n()
    if abs(error) >= 0.0001:
        print (n, "\t", trapezoidal, "\t", error)
    if abs(error) <= 0.0001:
        print(n, "nodes")

print("by the table above we need n=41 for s=1 to get error less than 1/100.")

trapezoidal = CompositeTrapezoidal(f, -1, 1, 41).n()
error = (trapezoidal - integral(f(x), x,-1,1)).n()
print("integral(f, x, -1, 1) ≈&quot;, trapezoidal, "by Composite Trapezoid Rule")
print()

#s=2---------------------------------------------------------------------------

print("For Standard Deviation, s, of 2, we have a=-2 and b=2")


def CompositeTrapezoidal(f, a, b, n):
    #compute the size of each subinterval
    h = (b-a)/n
    #here we add all the terms
    sum = 1/2*f(a)
    #add the values for all intermediate nodes
    for i in [1,2,..,n-1] :
        sum += f(a+i*h)
    #add the last node
    sum +=1/2*f(b)
    #multiply by h to complete the computation
    sum = sum*h
    #return the sum
    return sum
#define the function we will be using
f(x) =(e^-((x^2)/2))/(sqrt(2*pi))
#print the top part of the table
print ("n", "\t","trapezoidal" , "\t\t", "error")
print ("--------------------------------------------------")
#call the function several times
for power in [21,22,..,55]:
    n = power
    trapezoidal = CompositeTrapezoidal(f, -2, 2, n).n()
    error = (trapezoidal - integral(f(x), x,-2,2)).n()
    if abs(error) >= 0.0001:
        print (n, "\t", trapezoidal, "\t", error)
    if abs(error) <= 0.0001:
        print(n, "nodes")
print("By the table above we need n=54 when s=2 to get error less than 1/100.")

trapezoidal = CompositeTrapezoidal(f, -2, 2, 54).n()
error = (trapezoidal - integral(f(x), x,-2,2)).n()
print("integral(f, x, -2, 2) ≈&quot;, trapezoidal, "by Composite Trapezoid Rule")
print()




#s=3 ---------------------------------------------------------------

print("For Standard Deviation, s, of 3, we have a=-3 and b=3")

def CompositeTrapezoidal(f, a, b, n):
    #compute the size of each subinterval
    h = (b-a)/n
    #here we add all the terms
    sum = 1/2*f(a)
    #add the values for all intermediate nodes
    for i in [1,2,..,n-1] :
        sum += f(a+i*h)
    #add the last node
    sum +=1/2*f(b)
    #multiply by h to complete the computation
    sum = sum*h
    #return the sum
    return sum
#define the function we will be using
f(x) =(e^-((x^2)/2))/(sqrt(2*pi))
#print the top part of the table
print ("n", "\t","trapezoidal" , "\t\t", "error")
print ("--------------------------------------------------")
#call the function several times
for power in [21,22,..,55]:
    n = power
    trapezoidal = CompositeTrapezoidal(f, -3, 3, n).n()
    error = (trapezoidal - integral(f(x), x,-3,3)).n()
    if abs(error) >= 0.0001:
        print (n, "\t", trapezoidal, "\t", error)
    if abs(error) <= 0.0001:
        print(n, "nodes")
print("By the table above we need n=54 when s=3 to get error less than 1/100.")

trapezoidal = CompositeTrapezoidal(f, -3, 3, 54).n()
error = (trapezoidal - integral(f(x), x,-3,3)).n()
print("integral(f, x, -3, 3) ≈&quot;, trapezoidal, "by Composite Trapezoid Rule")

标签： nodesnormal-distributionsagestandard-deviationnumerical-integration