数值分析实验之数值积分法（Python 代码）

详细实验指导见上一篇，此处只写内容啦

实验内容

     选择 y=arctan(x) 在0-1上的积分

• 复化simpson算法

from sympy import *
import math

def func(x):
    return math.atan(x)
def sum_fun_xk(xk, func):
    return sum([func(each) for each in xk])

def integral(a, b, n, func):
    h = (b - a)/float(n)
    xk = [a + i*h for i in range(1, n)]
    return h/2 * (func(a) + 2 * sum_fun_xk(xk, func) + func(b))

if __name__ == "__main__":
     
    a, b = 0, 1
    n = 50
    y=integral(a, b, n, func)
    print ("复化simpson公式：",y)

运行结果：

    将区间分成50份所得结果：

      

    将区间分成100份所得结果：

     

• 复化梯形算法

import numpy 
from scipy import integrate
import matplotlib.pyplot as plt
import math

#x [0,2pi]
def f(x):
    return math.atan(x)

def T(a,b,n=50):
    h = (b - a) / n
    temp = 0
    for i in range(1,n):
        x = a + i * h
        temp += 2 * f(x)
    return (b - a) / (2 * n) * (f(a) + temp + f(b))

#50 100 200 500 1000
def S(a,b,n=50):
    h = (b - a) / n
    temp1 = 0
    temp2 = 0
    for i in range(1,n):
        xk1 = a + h * i
        xk2 = a + h * (i + 1)
        xk12 = (xk1 + xk2) / 2
        temp1 += f(xk1)
        temp2 += f(xk12)
    temp2 += f((a + a + h) / 2)
    return (b - a) / (6 * n) * (f(a) + 4 * temp2 + 2 * temp1 + f(b))

if __name__ == '__main__':
    n =1000 #50 100 200 500 1000
    a = 0
    b = 1

    print ('Truth-value:',integrate.quad(f,a,b)[0])
    print ('T-Estimated-value:',T(a,b,n))
    print ('S-Estimated-value:',S(a,b,n))

运行结果：

      

   

体会

  首先同一方法的比较，我们将区间分成了50份和100份进行计算，得出结果：将区间分成100份的计算结果精度较高。其次，不同方法之间的比较，我们选用的复化梯形算法和复化辛普森算法对该积分进行计算，发现在将区间分成50份时，所得精度相差不大。总而言之，通过本次实验，我对数值分析中相关的知识更加熟悉。其次，对各种编程软件的使用也进一步了解，虽然java是上学期学的，但是通过本次课进行熟悉，收获也颇多。