【Julia1.5.2】龙贝格积分

流星雨

参考数目：
【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.一般地，如果积分区间[a,b]太大或者f(x)变化较为剧烈，则可将[a,b]分段应用龙贝格积分。

Romberg.jl (2.6 KB, 下载次数: 1)

2020-12-26 19:48 上传

点击文件名下载附件

说明.docx (122.7 KB, 下载次数: 1)

2020-12-26 19:48 上传

点击文件名下载附件

watermelon

我也来个Py3.8.5的：

[Python] 纯文本查看 复制代码

# -*- 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[i+1]-self.x[i])*(self.y[i+1]+self.y[i]) / 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[1] + (T[1]-T[0])/3)
        S.append(T[2] + (T[2]-T[1])/3)
        S.append(T[3] + (T[3]-T[2])/3)
        
        C = []
        C.append(S[1] + (S[1]-S[0])/15)
        C.append(S[2] + (S[2]-S[1])/15)
        
        R = []
        print(T)
        print(self.N)
        R.append(C[1]+(C[1]-C[0])/63)
        # reserve the memory
        for i in range(self.__MAX_STEP):
            self.N = 2*self.N-1
            T[0:3] = T[1:4]
            T[3] = self.TrapeIntegral()
            S[0:2] = S[1:3]
            S[2] = T[3] + (T[3]-T[2])/3
            C[0] = C[1]
            C[1] = S[2] + (S[2]-S[1])/15
            R.append(C[1] + (C[1]-C[0])/63)
            if np.fabs(R[i+1] - R[i]) <= self.__ERR:
                print("Iteration has converged...\n")
                break
            print(T[3])
            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)