【Julia1.5.2】龙贝格积分
本帖最后由 流星雨 于 2020-12-26 19:49 编辑参考数目：
【1】《数值分析》(李乃成、梅立泉编著，科学出版社，ISBN：978-7-03-032192-3)
考虑到0xAA55大佬的提议，这里尽量避开繁杂的推导。具体亦可参考“https://blog.csdn.net/m0_37816922/article/details/103475445”
本程序暂时只能求解显函数积分。
需要用到的几个公式、定理和假设：
1.拉格朗日插值公式
2.广义积分中值定理
3.积分步长等距假设
4.在积分区间内高阶导数连续且变化不大
程序结构：
采用Julia语言(1.5.2)编制脚本Romberg.jl，主要功能为求解已知解析式的龙贝格数值积分。程序包含三个子函数：
1.fx(x)：定义函数解析式；
2.Σfx(a,b,k)：求函数的累加值，其中a为区间下界、b为区间上界，k为当前迭代次数；
3.ROMBERG_G(ε,fx)：龙贝格积分算法，其中ε为迭代偏差，fx为函数fx(x)；
用户只需更改函数fx(x)中的解析式和ROMBERG_G(ε,fx)中的迭代偏差ε。
使用说明:
1.安装Julia编译环境JuliaPro-1.5.2-1。可在URL：https://juliacomputing.com/products/juliapro/，下载最新版的JuliaPro-1.5.3-1；
2.安装完毕后用Juliapro打开文件Romberg.jl，点击▶即可运行。UI右侧Workspace窗口可实时显示各变量值；
3.其中answer为积分值、h为积分步长、number为迭代次数。用户可以更改a，b来调整积分区间，也可更改被积函数fx(x)。
附上参考书【1】给出的使用建议：
1.截断误差在步长h的8次方。但仍然要注意f^8 (η)的性态，若要用R_2n=R_n+(R_2n-R_n)/(4^4-1)，需保证f^10 (η)不要太大；
2.一般地，如果积分区间太大或者f(x)变化较为剧烈，则可将分段应用龙贝格积分。
我也来个Py3.8.5的：
# -*- coding: utf-8 -*-
"""
Created on Sat Dec 26 19:42:39 2020
Copyright@ Xi'an Jiaotong University - NuTHeL
@author: Zhou Xingguang
"""
import numpy as np
class Romberg(object):
def __init__(self, begin, end, N, func):
# critera condition
self.__MAX_STEP = 100
self.__ERR = 1E-12
self.begin = begin
self.end = end
self.N = N
self.func = func
# private method
def __GetPoint(self):
self.x = np.linspace(self.begin, self.end, self.N)
self.y = self.func(self.x)
def TrapeIntegral(self):
self.__GetPoint()
result = 0.0
for i in range(self.N-1):
result += (self.x-self.x)*(self.y+self.y) / 2.0
return result
def RombergIntegral(self):
R = 0
T = []
''' Initial Romberg '''
T.append(self.TrapeIntegral())
self.N = 2*self.N-1
T.append(self.TrapeIntegral())
self.N = 2*self.N-1
T.append(self.TrapeIntegral())
self.N = 2*self.N-1
T.append(self.TrapeIntegral())
S = []
S.append(T + (T-T)/3)
S.append(T + (T-T)/3)
S.append(T + (T-T)/3)
C = []
C.append(S + (S-S)/15)
C.append(S + (S-S)/15)
R = []
print(T)
print(self.N)
R.append(C+(C-C)/63)
# reserve the memory
for i in range(self.__MAX_STEP):
self.N = 2*self.N-1
T = T
T = self.TrapeIntegral()
S = S
S = T + (T-T)/3
C = C
C = S + (S-S)/15
R.append(C + (C-C)/63)
if np.fabs(R - R) <= self.__ERR:
print("Iteration has converged...\n")
break
print(T)
print(self.N)
return R
def func(x):
return np.sin(x)/x
if __name__ == '__main__':
r = Romberg(1e-16, 1, 2, func)
result = r.RombergIntegral()
print(result)
页:
[1]