Romberg 積分的python 程式

2021-08-31 09:39:39 字數 2251 閱讀 1038

程式作用:用romberg演算法計算一重積分

語言:python

具體演算法見鏈結

語言:python

程式介紹

本程式包含三個函式

t_2n:計算t(h/2)

romberg:積分主程式

fun(x):函式形式

另外還有乙個主函式

romberg函式的演算法:

該函式中採用了兩個陣列來分別儲存相鄰兩行的值,這樣做的優點是演算法簡單清晰,計算速度快,缺點是占用儲存空間。該函式使用了兩重迴圈,第一重迴圈計算第一列t0,第二重迴圈計算每一行的t(m+1), 最後通過比較對角線上數的差值以及相鄰兩個數的差值來控制迴圈結束,為防止出現死迴圈設定了最大迭代次數k

直接上**

import math

import numpy as np

from numpy import

*def

t_2n

(a, b, n, t_n)

:#計算t0(h/2)的函式, a 積分下限,b 積分上限, n是區間等分數

if n<1:

print

('n should larger than 1'

) h =

(b - a)

/n #步長

sum_f =0.

#初始化,中間變數

for k in

range(0

, n)

: sum_f = sum_f + fun(a +

(k +

0.5)

*h) t_2n = t_n/2.

+ sum_f*h/2.

return t_2n

defromberg

(a, b, err_min)

: kmax =

6 tm = zeros(kmax,dtype =

float

)# 第m行所有的元素

tm1 = zeros(kmax,dtype =

float

)#第m+1行所有的元素

tm[0]

=0.5

*(b-a)

*(fun(a)

+ fun(b)

)# 初始值

print

(tm)

err =1.

k =0 np.set_printoptions(precision =9)

while

(err>err_min and k

:#控制迴圈次數

n =2**k # n是區間等分數

m =1 tm1[0]

= t_2n(a, b, n, tm[0]

)while

(err>err_min and m <=

(k+1))

:#控制迴圈次數

tm1[m]

= tm1[m-1]

+(tm1[m-1]

-(tm[m-1]

))/(

4.**m-1)

result = tm1[m]

err1 =

abs(tm1[m]

-tm[m-1]

) err2 =

abs(tm1[m]

-tm1[m-1]

) err =

min(err1,err2)

m = m+

1 tm = np.copy(tm1)

k = k+

1print

(tm)

return result

deffun

(x):

#被積函式, x為自變數

if x <=0:

print

('err'

)return

99999

else

: f = math.log(x)

return f

f1 = romberg(1,

2,1.e-12

)#主程式(分別輸入上下限,誤差)

print

('result = '

,f1)

最終結果:

數值積分 龍貝格 Romberg 積分

數值積分在工程上是個比較有用的數學工具。在工程上有很多數學問題,看似簡單,計算所用的數學公式不算複雜,但是求解起來卻很困難,很難獲得解析解的公式,這個時候就需要用到數值求解的辦法來獲取滿足工程需要的近似數值結果。舉個簡單的例子,求乙個變截面圓柱 或者叫圓台 的等效截面慣性矩,計算公式為ie h 3 ...

用C 寫的龍貝格(Romberg)積分法

很久之前了,群上有人問有沒有用c 寫的積分函式,因為我自己以前做過乙個積分計算器,就跟他說搜尋一下romberg積分法吧,但是後來他說網上寫的都是用c 寫的方法,沒有c 的,問我有沒有,但是我之前是在家裡做的,學校裡沒有副本,可惜。今天翻一下郵箱,發現居然是我忘記貼標籤,導致當時沒有找到,感覺有點對...

積分學 重積分與曲線積分曲面積分的理解

積分作為高等數學的核心部分,主要含蓋了一重積分,二重積分,三重積分,第一型曲線積分,第二型曲線積分,第一型曲面積分,第二型曲面積分。微積分學在研究中作為必不可少的工具,熟練掌握一些計算方法和重要公式比如是最基本的了。下面是我的一些總結 一重積分,主要精力就要研究不定積分和定積分了。不定積分的求解是後...