技术宅的结界

 找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 115|回复: 1
收起左侧

【Julia1.5.2】龙贝格积分

[复制链接]

4

主题

7

帖子

360

积分

用户组: 中·技术宅

UID
6285
精华
3
威望
42 点
宅币
203 个
贡献
51 次
宅之契约
0 份
在线时间
4 小时
注册时间
2020-10-2
发表于 2020-12-26 19:37:55 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有帐号?立即注册→加入我们

x
本帖最后由 流星雨 于 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.在积分区间内高阶导数连续且变化不大

01.png 02.png

程序结构:
采用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)中的迭代偏差ε。

03.png

使用说明:
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)
回复

使用道具 举报

28

主题

308

帖子

1811

积分

用户组: 上·技术宅

UID
3808
精华
10
威望
99 点
宅币
1100 个
贡献
155 次
宅之契约
0 份
在线时间
332 小时
注册时间
2018-5-6
发表于 2020-12-26 23:13:14 | 显示全部楼层
我也来个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)
Passion Coding!

本版积分规则

QQ|申请友链||Archiver|手机版|小黑屋|技术宅的结界 ( 滇ICP备16008837号 )|网站地图  

GMT+8, 2021-1-15 22:36 , Processed in 0.092591 second(s), 31 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表