复化求积公式

来源:互联网 发布:程序员的头像 编辑:程序博客网 时间:2024/04/30 05:01

一、简介

复化求积公式(composite integration rule )一类重要的求积公式,指将求积区间分为n个子区间,对每个子区间应用同一求积公式,所得到的复合数值积分公式。

详解

二、实现

# -*- coding: utf-8 -*-"""Created on Mon Dec 19 10:21:15 2016    复化梯形、复化Simpson@author: Administrator"""from numpy import *from scipy import integrateimport matplotlib.pyplot as plt#x [0,2pi]def f(x):    return e**(3*x) * cos(pi * 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 1000def 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 = 2 * pi    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)    x = linspace(0,2*pi,100)    y = f(x)    fig = plt.figure(figsize=(8,4))    ax = fig.add_subplot(111)    ax.set_xlabel('x')    ax.set_ylabel('y')    ax.plot(x,y,color='r',linestyle='-',label='f(x)')    ax.legend(loc='upper left')    fig.show()    fig.savefig('a.png')

原函数

0 0
原创粉丝点击