计算科学 - 改变积分变量以避免无穷大 - 吾爱随笔录

改变积分变量以避免无穷大

计算科学数字龙格库塔一体化

2021-11-30 05:12:52

我想在 Fortran 中编写代码以数值方式求解此积分：最好的方法是什么？

\int_{1095}^{\infty} \frac{d x}{x \sqrt{(x + 644.153) (4.17 \cdot 10^{- 5} x + 0.145)}}

$\int_{1095}^\infty \frac{dx}{x\sqrt{(x+644.153)(4.17 \cdot 10^{-5} x+0.145)}}$

我尝试了 Monte-Carlo 和 RK4，但我不得不输入一个大数字而不是无穷大。

为了避免无穷大，我将变量更改为。那么积分变为： $x$ $1/t$

\int_{0}^{1 / 1095} \frac{d t}{(t \sqrt{(1 / t + 644.153) (4.17 \cdot 10^{- 5} / t + 0.145)}}

$\int_0^{1/1095}\frac{dt}{(t\sqrt{(1/t+644.153)(4.17 \cdot 10^{-5}/t+0.145)}}$

但这一次的问题是，积分的下限是分母的根之一。你能建议另一种改变变量的方法来避免这个问题吗？

3个回答

你快到了，只要把放在平方根下，它就会变成 $t$

$\displaystyle\int_0^{1/1095} {\frac{dt}{\sqrt{(1+644.153t)(4.17 \cdot 10^{-5} + 0.145 t)}}},$

这消除了 $t=0$

有一个封闭形式的解决方案。为了 $a,b,c,d,x_0 > 0$ ,

\int_{x_{0}}^{\infty} \frac{d x}{x \sqrt{(a + b x) (c + d x)}} = \frac{\log (2 \sqrt{b d (a + b / x_{0}) (c + d / x_{0})} + a d + b c + 2 b d / x_{0}) - \log (2 \sqrt{b d a c} + a d + b c)}{\sqrt{b d}}

$\int_{x_0}^\infty \frac{dx}{x\sqrt{(a+bx)(c+dx)}} = \\ \frac{\log(2 \sqrt{b d (a + b/x_0) (c + d/x_0)} + a d + b c + 2 b d/x_0) - \log(2 \sqrt{b d a c} + a d + b c)}{\sqrt{b d}}$

a <- 1
b <- 644.153
c <- 4.17e-5
d <- 0.145
x <- 1/1095
(log(2*sqrt(b*d*(a+b*x)*(c+d*x)) + a*d + b*c + 2*b*d*x) - 
    log(2*sqrt(a*b*c*d) + a*d + b*c))/sqrt(b*d) 
# 0.08334367

使用 $c = 4.15e^{-5}$ 作为@nicoguaro，我们发现与他的结果相同：

c <- 4.15e-5
(log(2*sqrt(b*d*(a+b*x)*(c+d*x)) + a*d + b*c + 2*b*d*x) - 
    log(2*sqrt(a*b*c*d) + a*d + b*c))/sqrt(b*d) 
# 0.08344436

我认为你不需要改变这个问题的变量（你自己）。Quadpack似乎很适合它。它使用Gauss-Kronrod正交。

我尝试过（在 Python 中），它似乎有效。

import numpy as np
from scipy.integrate import quad

fun = lambda x: 1/(x*np.sqrt((x + 644.153)*(4.15e-5*x + 0.145)))
inte, err = quad(fun, 1095, np.inf)

这给出0.08344435887373103了错误的结果2.7989753696121526e-10。

其它你可能感兴趣的问题

上一篇实施[ X, ⋅ ][X,⋅]作为一个矩阵，其中是一个矩阵n2×n2n2×n2XXn × nn×n 下一篇为什么这个 3D ising 自旋系统的磁化表现出突然的行为