【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的：

1. # -*- coding: utf-8 -*-
   
2. """
   
3. Created on Sat Dec 26 19:42:39 2020
   
4. Copyright@ Xi'an Jiaotong University - NuTHeL
   
5. @author: Zhou Xingguang
   
6. """
   
7. 
   
8. import numpy as np
   
9. 
   
10. class Romberg(object):
    
11.     def __init__(self, begin, end, N, func):
    
12.         # critera condition
    
13.         self.__MAX_STEP = 100
    
14.         self.__ERR = 1E-12
    
15.         
    
16.         self.begin = begin
    
17.         self.end = end
    
18.         self.N = N
    
19.         self.func = func
    
20.    
    
21.     # private method
    
22.     def __GetPoint(self):
    
23.         self.x = np.linspace(self.begin, self.end, self.N)
    
24.         self.y = self.func(self.x)
    
25.    
    
26.     def TrapeIntegral(self):
    
27.         self.__GetPoint()
    
28.         result = 0.0
    
29.         for i in range(self.N-1):
    
30.             result += (self.x[i+1]-self.x[i])*(self.y[i+1]+self.y[i]) / 2.0
    
31.         return result
    
32.    
    
33.     def RombergIntegral(self):
    
34.         R = 0
    
35.         T = []
    
36.         ''' Initial Romberg '''
    
37.         T.append(self.TrapeIntegral())
    
38.         self.N = 2*self.N-1
    
39.         T.append(self.TrapeIntegral())
    
40.         self.N = 2*self.N-1
    
41.         T.append(self.TrapeIntegral())
    
42.         self.N = 2*self.N-1
    
43.         T.append(self.TrapeIntegral())
    
44.         
    
45.         S = []
    
46.         S.append(T[1] + (T[1]-T[0])/3)
    
47.         S.append(T[2] + (T[2]-T[1])/3)
    
48.         S.append(T[3] + (T[3]-T[2])/3)
    
49.         
    
50.         C = []
    
51.         C.append(S[1] + (S[1]-S[0])/15)
    
52.         C.append(S[2] + (S[2]-S[1])/15)
    
53.         
    
54.         R = []
    
55.         print(T)
    
56.         print(self.N)
    
57.         R.append(C[1]+(C[1]-C[0])/63)
    
58.         # reserve the memory
    
59.         for i in range(self.__MAX_STEP):
    
60.             self.N = 2*self.N-1
    
61.             T[0:3] = T[1:4]
    
62.             T[3] = self.TrapeIntegral()
    
63.             S[0:2] = S[1:3]
    
64.             S[2] = T[3] + (T[3]-T[2])/3
    
65.             C[0] = C[1]
    
66.             C[1] = S[2] + (S[2]-S[1])/15
    
67.             R.append(C[1] + (C[1]-C[0])/63)
    
68.             if np.fabs(R[i+1] - R[i]) <= self.__ERR:
    
69.                 print("Iteration has converged...\n")
    
70.                 break
    
71.             print(T[3])
    
72.             print(self.N)
    
73.         return R
    
74.    
    
75. 
    
76. def func(x):
    
77.     return np.sin(x)/x
    
78. 
    
79. 
    
80. if __name__ == '__main__':
    
81.     r = Romberg(1e-16, 1, 2, func)
    
82.     result = r.RombergIntegral()
    
83.     print(result)
    

复制代码