解決不正確的積分不近似

我在解決python中的這個積分問題。集成的功能沒有在集成的邊界上定義。我發現了幾個類似於這個問題的更多問題，但所有問題都是非常具體的回答。我不想近似積分太多，如果可能的話，根本就不是，因爲我在做這個積分首先是爲了避免近似。有什麼辦法解決這個積分？解決不正確的積分不近似

import numpy as np 
from pylab import * 
import scipy 
from math import * 
from scipy import integrate 

m_Earth_air = (28.0134*0.78084)+(31.9988*0.209476)+(39.948*0.00934)+(44.00995*0.000314)+(20.183*0.00001818)+(4.0026*0.00000524)+(83.80*0.00000114)+(131.30*0.000000087)+(16.04303*0.000002)+(2.01594*0.0000005) 
Tb0 = 288.15 
Lb0 = -6.5 
Hb0 = 0.0 
def Tm_0(z): 
    return Tb0+Lb0*(z-Hb0) 
k = 1.38*10**-19 #cm^2.kg/s^2.K #Boltzmann cst 
mp = 1.67262177*10**-27 #kg 
Rad= 637100000.0 #radius planet #cm 
g0 = 980.665 #cm/s^2 
def g(z): 
    return (g0*((Rad/(Rad+z))**2.0)) 
def scale_height0(z): 
    return k*Tm_0(z*10**-5)/(m_Earth_air*mp*g(z)) 



def functionz(z,zvar): 
    return np.exp(-zvar/scale_height0(z))*((Rad+zvar)/(Rad+z))/((np.sqrt(((Rad+zvar)/(Rad+z))**2.0-1.0))) 

def chapman0(z): 
    return (1.0/(scale_height0(z)))*((integrate.quad(lambda zvar: functionz(z,zvar), z, np.inf))[0]) 

print chapman0(1000000) 
print chapman0(5000000)

第一塊變量和定義都很好。問題在於「functionz（z，zvar）」及其整合。非常感謝任何幫助！

來源

2015-06-11 JadeChee

您是否期待SciPy通過分析解決這個問題？它不這樣做; scipy.integrate例程都計算數值近似值。 – user2357112

'm_Earth_air'的價值是什麼？ – fjarri

它通常通過重新縮放變量來消除可能的數值不穩定性。在你的情況下，zvar從1e6開始，這可能是由於quad()中的某些實現細節引起的問題。如果您縮放它作爲y = zvar/z，使整合從1開始似乎收斂相當不錯z = 1e6：

def functiony(z, y): 
    return np.exp(-y*z/scale_height0(z))*(Rad+y*z)/(Rad+z)/np.sqrt(((Rad+y*z)/(Rad+z))**2.0-1.0) 

def chapman0y(z): 
    return (1.0/(scale_height0(z)))*((integrate.quad(lambda y: functiony(z,y), 1, np.inf))[0]) 

>>> print(chapman0y(1000000)) 

1.6217257661844094e-06

（我設置m_Earth_air = 28.8e-3 - 這個常量在代碼中失蹤，我以爲這是摩爾質量（編輯）kg/mole）。

至於z = 5e6,scale_height0(z)是負數，它在指數下給出了一個巨大的正值，使得積分在無窮大上發散。

來源

2015-06-11 05:47:50 fjarri

儘管它現在發佈了一個值，但問題仍然是當y = 1.0時函數（z，y）是無限的，當y = inf時是nan。正因爲如此，積分的價值比它應該小得多。 Ps。你對scale_height0（z）是正確的，實際上我對z的不同範圍有不同的標度高度函數，而scale_height0（z）只能用於高達11km（我的單位是11e5）。 Ps2。很好的猜測m_Earth_air！它是28.8千克/千摩爾，平均空氣分子量:) – JadeChee

'quad（）'在積分收斂時對無窮無問題。你可以嘗試將'exp（-x）/ sqrt（x）'從'0'集成到'inf'，你將得到'sqrt（pi）'。你的積分也沒有任何不妥之處。你怎麼知道這個值比它應該小？ Wolframalpha集成它沒有問題，並給出了相同的答案。 – fjarri

除非你可以解析積分，否則沒有辦法解決它，如果沒有一個近似的邊界。這不是一個Python問題，而是一般的微積分問題，因此數學課爲什麼要用這麼大的努力向你展示數值近似值。

如果你不希望它差異太大，選擇一個小的epsilon與快速收斂的方法。

編輯 - 清晰度上最後一條語句：

的Epsilon - ɛ - 指的是通過集成 - 三角洲的邊界步長的x記住數字的近似方法所有切片的積分成條子並將它們添加回將其視爲每個條子的寬度，條子越小，近似值越好。您可以在數字包中指定它們。

快速收斂的方法意味着該方法快速逼近積分的真值，並且對於每個條子，近似誤差都很小。例如，黎曼和是一種天真的方法，它假定每個條子都是矩形，而梯形將條子的開始和結束與一條線連接以形成梯形。在這兩個中，梯形通常會更快地收斂，因爲它試圖解釋形狀內的變化。（因爲大多數函數都有更好的猜測，所以通常都不使用）

這兩個變量都會改變計算的計算開銷。通常epsilon是最昂貴的改變，因此爲什麼選擇一個好的近似方法很重要（對於同一個epsilon，有些差異可能會有一個數量級的誤差）。

所有這些都取決於您的計算可以容忍多少錯誤。

來源

2015-06-11 04:18:52 John

當你說「用快速收斂的方法選擇一個小epsilon」時，你在想什麼方法？ – JadeChee

我有類似的問題，發現SciPy quad需要你指定另一個參數，epsabs=1e-1000,limit=1000（stepsize limit），epsrel=1e1適用於我試過的所有東西。即在這種情況下：

def chapman0(z):  
    return (1.0/(scale_height0(z)))*((integrate.quad(lambda zvar: functionz(z,zvar), z, np.inf, limit=1000, epsabs=1e-1000, epsrel=1e1))[0])[0]) 
#results: 
0.48529410529321887 
-1.276589093231806e+21

似乎是一個絕對高容錯性，但對於不快速收斂，似乎解決這個問題積分。只是張貼爲其他類似問題的人，因爲這篇文章是相當過時的。其他軟件包中的算法收斂速度更快，但在SciPy中找不到。結果基於發佈的代碼（不是選定的答案）。

來源

2016-06-08 00:58:10 Matt

解決不正確的積分不近似

回答

相關問題