Numerical integration in Python with unknown constant

Question

I’d like to solve the below equation for the unknown $T$:

$$\int_0^\infty \frac{x^2}{\exp\left(\frac{x}{T}\right)-1}\kappa_x \mathrm{d}x = C,$$

where $C$ is a known constant and $\kappa_x$ is some function of $x$, although the actual function is unknown – all I have as an array of values of $\kappa$ over different values of $x$. I do not want to have to fit this it to a function first. I was hoping to use one of the SciPy’s numerical integration functions such as integrate.simps.

I’m confused on how to solve this; since $T$ is inside an exponential with $x$ dependence I cannot simply pull it out of the integral. It would be great to be able to pass the unknown constant into my numerical integral and receive an answer in terms of the constant. But as far as I know, in Python, I cannot pass an unknown. How can I accomplish this?

Edit: Here is an example of what my $\kappa_x$ array could look like:

k_arr = [1.1e-3, 8.8e-2, ..., 3.7e3]

(i.e., floating point numbers generally ranging from ~ 10^-3 to 10^3. But the length of my real array is about 3000).

To provide some more explanation, the reason I subscript it as $\kappa_x$ is because each of the values is associated with some $x$ value (in my case, $x$ is an energy or frequency). Previously I used various energies to arrive at $\kappa$ values, which are opacities in units of cm^2 g^-1. The integral in question relates to the power per unit mass emitted by dust in a galaxy, where I omitted constants (except for the temperature $T$). It is equal to some known number, so that I can hopefully solve for the temperature of the dust.

Answer 1

What you want seems inherently impossible, and that’s not due to restrictions of Python.

The only way we can arrive at a situation where we only need to apply a single quadrature is to get analytically get rid of all dependencies of $T$ in the integral. To this end, the best we can do is to apply the substitution $x=Ty$ to your integral: $$ I(T) = \int_0^∞ \frac{x^2 κ_x}{\exp\left(\frac{x}{T}\right)-1} \mathrm{d}x = \int_0^∞ \frac{T^2 y^2 κ_{Ty}}{\exp(y)-1} T \mathrm{d}y = T^3\int_0^∞ \frac{y^2 κ_{Ty}}{\exp(y)-1} \mathrm{d}y. $$ Now, we got rid of all contributions of $T$ to the integral except as an argument of $κ$. We can only improve this further by applying special knowledge about $κ$, but since $κ$ is empirical, you probably do not know anything that could help you here.

From another perspective, if you apply a method similar to the midpoint method to the values of $κ$ that you actually have ($κ_{x_1},…,κ_{x_n}$), the result $I(T)$ would be some weighted sum of these values: $$ I(T) ≈ \sum_{i=1}^n a_i(T) κ_{x_i}, $$ with the weights for the being: $$ a_i(T) = \frac{x_i^2}{\exp\left(\frac{x_i}{T}\right)-1} b_i,$$ where $b_i$ only depends on the spacing of your $x_i$. Obviously, you cannot just factor out the $T$ dependency here. For other quadrature methods, the weights may become more complicated, but there is no reason to expect this solves the above problem.

Answer 2

It seems to me that this question is about how you handle $\kappa$. I know you don't want to fit a function to it, but I don't think you can get around sampling it at positions other than the positions in the array. You may also need some kind of asymptotic approximation as $x \rightarrow \infty$ so you can extrapolate in that range.

Having said that, I have an idea about how to handle the dependency on $T$ that may be easier.

I don't know anything about SciPy, so you will need to do some algebra.

The key observation that you need is: If $T > 0$ and $k$ is a positive integer then:

$$\int_{0}^{\infty} \frac{x^k}{e^{x/T} -1}\, dx = k!\,T^{k+1} \zeta(k+1)$$

Where $\zeta(z)$ is the Riemann zeta function. (Exercise: Prove this.)

The values of the zeta function at positive integers are well-known. This gives you a way to calculate:

$$\int_{0}^{\infty} p(x) \frac{x^2}{e^{x/T} -1}\, dx$$

for any polynomial $p$. So if $\kappa$ were a polynomial, you could calculate this more or less directly to obtain a polynomial in $T$, and then solve for $T$.

But if $\kappa$ is conceptually defined in the range $[0,\infty)$, chances are there isn't a good polynomial approximation for it that behaves sensibly over the whole range.

So my suggestion is a variant of Wrzlprmft's, but instead of using the midpoint method, use Gaussian quadrature.

If you find the Gauss weights and nodes (which will be functions of $T$) for the weight function $W_T(x) = \frac{x^2}{e^{x/T} - 1}$ over the domain $[0,\infty)$, this will give you an approximation to the integral expressed as a finite set of samples of $\kappa$.

This may give you a well-behaved equation that you can solve for $T$.

Answer 3

Lets assume that the given expression has a non-infinite integral you could actually write down.

The question you should be asking -might- be more like: "To what precision do I need to calculate?". If you have $$K_x(x)$$ only as an array of values anyway, you might only ever achieve a certain precision. Is an error of 1e-16 okay for your application? Do you need even more precision?

You might approach the problem from that perspective by starting to evaluate the expression at equally spaced points x (e.g. the spacing of the array of K_x you have), and apply the weights like @Wrzlprmft wrote. You should be able to estimate how big your error is if you stop the summation at a certain point. The exponential function will force the profile towards zero quite rapidly, so when you are at a high value of x_i of your integration points, you probably already reach the smallest value storable in double precision. Depending on the problem you want to solve that might already be enough. Could you elaborate what the underlying problem is?

Answer 4

This could be done purely numerically, if you take the integral as part of your non-linear equation.

First of all, I created a interpolating function for your $\kappa_x$ function. It's just a dummy linear interpolation between data points allowing for end-point constant extrapolation (see interp1d). I leave it up to you to modify this according to your problem and your taste (more fancy interpolation scheme).

Then I define the left-hand-side of the non-linear equation $f(T) = 0$ as $$f(T) = \int_{0}^{\infty} \frac{x^{2}}{\exp(\frac{x}{T})-1}\kappa_{x}\, dx -C.$$ So I need to define a function that takes only takes $T$ as an argument and encapsulate the integral in this function. In Python, this could be something like the code below. Note that I integrate from the machine epsilon on the left boundary to avoid the singularity at $x = 0$.

Make sure you read up on the docs for quad and fsolve so that you understand their options and limitations. Depending on your $\kappa_{x}$ function, your integrand could be peaked.

import numpy as np
from scipy.integrate import quad
from scipy.optimize import fsolve
from scipy.interpolate import interp1d

kappa_x = np.linspace(0, 10, 11)
kappa_k = np.array([1, 2, 3, 4, 5, 5, 4, 3, 2, 1, 1])

kappa_fun = interp1d(kappa_x, kappa_k, bounds_error = False, fill_value = 1.0)

def RHS(T):
    C = 1000
    def integrand(x):
        return np.power(x, 2) / (np.exp(x / T) - 1.0) * kappa_fun(x)
    res, error = quad(integrand, np.finfo(float).eps, np.inf)
    return res - C

T_sol = fsolve(RHS, 100)