python辛普森积分_Python龙贝格法求积分实例

令Tn为将[a,b]划分n等分的复合梯形求积公式，h =(b-a)/n为小区间的长度。h/2类似于梯形公式中的(b-a)/2

注意：这里的k+1是下标

通过研究我们发现：T2n与Tn之间存在一些递推关系。

注意：这里的k+1/2是下标。并且其中的h/2是中的h是Tn(n等分中的h = (b-a)/n))

于是乎，我们可以一次推出T1,T2,T4,T8…T2n序列

引出这些之后，才是我们的主题：龙贝格求积公式

龙贝格求积公式的实质是用T2n序列构造，S2n序列，

再用S2n序列构造C2n序列

最后用C2n序列构造R2n序列。

编程实现，理解下面的几个公式即可。

python编程代码如下：

# coding=UTF-8

# Author:winyn

'''

给定一个函数，如：f(x)= x^(3/2)，和积分上下限a,b，用机械求积Romberg公式求积分。

'''

import numpy as np

def func(x):

return x**(3/2)

class Romberg:

def __init__(self, integ_dowlimit, integ_uplimit):

'''

初始化积分上限integ_uplimit和积分下限integ_dowlimit

输入一个函数，输出函数在积分上下限的积分

'''

self.integ_uplimit = integ_uplimit

self.integ_dowlimit = integ_dowlimit

def calc(self):

'''

计算Richardson外推算法的四个序列

'''

t_seq1 = np.zeros(5, 'f')

s_seq2 = np.zeros(4, 'f')

c_seq3 = np.zeros(3, 'f')

r_seq4 = np.zeros(2, 'f')

# 循环生成hm间距序列

hm = [(self.integ_uplimit - self.integ_dowlimit) / (2 ** i) for i in range(0,5)]

print(hm)

# 循环生成t_seq1

fa = func(self.integ_dowlimit)

fb = func(self.integ_uplimit)

t0 = (1 / 2) * (self.integ_uplimit - self.integ_dowlimit) * (fa+fb)

t_seq1[0] = t0

for i in range(1, 5):

sum = 0

# 多出来的点的累加和

for each in range(1, 2**i,2):

sum =sum + hm[i]*func( self.integ_dowlimit+each * hm[i])#计算两项值

temp1 = 1 / 2 * t_seq1[i - 1]

temp2 =sum

temp = temp1 + temp2

# 求t_seql的1-4位

t_seq1[i] = temp

print('T序列：'+ str(list(t_seq1)))

# 循环生成s_seq2

s_seq2 = [round((4 * t_seq1[i + 1] - t_seq1[i]) / 3,6) for i in range(0, 4)]

print('S序列：' + str(list(s_seq2)))

# 循环生成c_seq3

c_seq3 = [round((4 ** 2 * s_seq2[i + 1] - s_seq2[i]) / (4 ** 2 - 1),6) for i in range(0, 3)]

print('C序列：' + str(list(c_seq3)))

# 循环生成r_seq4

r_seq4 = [round((4 ** 3 * c_seq3[i + 1] - c_seq3[i]) / (4 ** 3 - 1),6) for i in range(0, 2)]

print('R序列：' + str(list(r_seq4)))

return 'end'

rom = Romberg(0, 1)

print(rom.calc())

以上这篇Python龙贝格法求积分实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持python博客。