Całka w przestrzeni log-log

Pracuję z funkcjami, które generalnie są znacznie płynniejsze i lepiej zachowują się w przestrzeni log-log --- więc tam wykonuję interpolację / ekstrapolację itp. I to działa bardzo dobrze. Czy istnieje sposób na zintegrowanie tych funkcji numerycznych w przestrzeni dziennika?

tzn. mam nadzieję, że użyję jakiejś prostej reguły trapezoidalnej, aby wykonać całkę skumulowaną (np. w pythonie, użyj scipy.integrate.cumtrapz), aby znaleźć st$F(r)$

$$F(r) = \int_0^r y(x) \, dx$$

Mam jednak nadzieję, że użyję wartości i zamiast i (jeśli to możliwe).$log(y)$$log(x)$$y$$x$

Możesz po prostu zmienić zmienne. Oprawa$a=log(x)$, $b(a)=log(y(x))$. Całka staje się

$F(r)=\int^{log(r)}_{-\infty} exp(a+b) da$

Musisz być trochę ostrożny, ponieważ integrujesz się z $-\infty$. To, co musisz dokładnie zrobić, będzie zależeć od czego$y(x)$ wygląda jak.

Nie używam pytona, ale jeśli dobrze rozumiem, to przez

$$F(r)=\int_0^ry(x)dx$$ myślisz coś takiego $$\mathbf{F}=\mathrm{integrate}(\mathbf{y},\mathbf{x})$$ gdzie $\mathbf{F}=[F_1,...,F_n]$ jest wektorem próbkującym całkę na siatce $\mathbf{x}$.

Jednak nie masz próbek $x$ i $y$, ale raczej masz próbki $\hat{x}=\log(x)$ i $\hat{y}=\log(y)$.

Oczywiście najprostsze podejście byłoby

$$\mathbf{F}=\mathrm{integrate}(\exp(\mathbf{\hat{y}}),\exp(\mathbf{\hat{x}})) ,$$ ale byłoby to podatne na błędy, ponieważ $y(x)$ mimo to nie jest gładka $\hat{y}(\hat{x})$ jest.

Teraz reguła trapezoidalna zasadniczo zakłada twój wkład$y(x)$jest fragmentarycznie liniowy. Zatem proste uogólnienie byłoby takie założenie$\hat{y}(\hat{x})$ jest fragmentarycznie liniowy.

W tym przypadku definiowanie $\Delta F_k=F_{k+1}-F_k$, ty masz

$$\Delta F_k=\int_{x_k}^{x_{k+1}}y(x)dx =\int_{\hat{x}_k}^{\hat{x}_{k+1}}e^{\hat{y}(\hat{x})}e^\hat{x}d\hat{x} =\int_{\hat{x}_k}^{\hat{x}_{k+1}}\tilde{y}(\hat{x})d\hat{x}$$

Następnie definiowanie $t=(\hat{x}-\hat{x}_k)/\Delta \hat{x}_k$, ty masz

$$\hat{y}_{k+t} \approx \hat{y}_k + t \Delta \hat{y}_k$$ i $\tilde{y}(t) \approx a e^{bt}$, z $a=e^{\hat{y}_k+\hat{x}_k}$ i $b=\Delta \hat{y}_k+\Delta \hat{x}_k$.

Tak więc całka staje się

$$\Delta F_k \approx a \Delta \hat{x} \int_0^1e^{bt}dt = a \Delta \hat{x} \frac{e^b-1}{b}$$

W Matlabie wyglądałoby to mniej więcej tak

dlogx=diff(logx); dlogy=diff(logy); k=1:length(logx)-1;  
b=dlogx+dlogy; a=exp(logx+logy);  
dF=a(k).*dlogx.*(exp(b)-1)./b;  
F=cumsum([0,dF]);

Mam nadzieję że to pomoże!

(Edycja: Moja odpowiedź jest zasadniczo taka sama, jak o wiele bardziej zwięzła odpowiedź, którą Damascus Steel podał podczas pisania. Jedyna różnica polega na tym, że próbowałem podać konkretne rozwiązanie dla przypadku, w którym „ $y(x)$„jest fragmentarycznie liniowy $\hat{y}(\hat{x})$ dyskretnie dyskretnie $\mathbf{\hat{x}}$ siatka, z $F(\hat{x}_1)=0$.)

Jeśli funkcja wygląda gładko na wykresie log-log, można interpolować za pomocą prawa mocy dla każdego przedziału (prawa mocy są oczywiście liniowe w log-log). Tak więc między punktami$(x_i, y_i)$ i $(x_{i+1}, y_{i+1})$ przy założeniu, że $y = C_i x^{n_i}$ w przedziale czasowym $i$, otrzymujesz $n_i = \ln(y_i/y_{i+1})/\ln(x_i/x_{i+1})$ and $C_i = \ln(y_i) - n_i\ln(x_i)$. The contribution to the integral from interval $i$ is then

$$\Delta F_i = \int_{x_i}^{x_{i+1}} C_i x^{n_i}dx = \left\{ \begin{array}{l@{\quad}l} \frac{C_i}{n_i + 1} (x_{i+1}^{n_i+1} - x_i^{n_i+1}), & n_i \neq -1 \\ C_i(\ln x_{i+1} - \ln x_i), & n_i=-1, \end{array} \right.$$ where you obviously need some tolerance for identifying the special case $n_i=-1$ in your implementation.

I think that there is a bit of confusion with change of variables in some of the previous answers as well as some errors. The integral of a log function is not the log of the integral. I think in general is difficult to write out the the integral of a function knowing the integral of its log. If anyone knows how to do that I would be interested.

In the meanwhile, @Stefan's solution above is the way to get around integrating a function in log-log space. The starting point is that the function you're dealing is linear in log-log space for small enough segments.

One can then write the equation of the line at the segment endpoints:

$$\log(y_1)=m_1 \log(x_1) + n_1$$ $$\log(y_2)=m_1 \log(x_1) + n_1$$

where $m_1$ is the slope of the line and $n_1$ is its y-intercept.

By subtracting the two, one can find:

$$m_1=\frac{\log(y_1) -log(y_2)}{\log(x_1)-log(x_2)}$$

i od podstawienia:

$$n_1=\log(y_1)-m_1 \log(x_1)$$

Jeśli w przestrzeni log-log równanie segmentu jest zbliżone do linii, to w normalnej (liniowej) przestrzeni równanie segmentu jest bliskie wykładniczemu:

$$y(x)\simeq x^{m} e^{n}$$

Jeśli mamy formułę analityczną dla tego segmentu, łatwo jest zintegrować:

$$\int_{x_1}^{x_2}y(x)dx = \frac{e^{n_1}}{m_1+1}(x_2^{m_1+1}-x_1^{m_1+1}), \,\, \text{for } \, m\neq -1$$

i

$$\int_{x_1}^{x_2}y(x)dx = e^{n_1} \log \frac{x_2}{x_1},\, \,\text{for } \, m = -1$$

To trochę przypomina oszustwo, ale jest to próbkowanie w przestrzeni log-logu, dzięki czemu możemy przybliżyć funkcję w przestrzeni liniowej do wykładniczej z parametrami pochodzącymi z przestrzeni log-log.

Rozwiązanie, którego używam, jest w zasadzie implementacją reguły trapezu i wykorzystuje scipy.misc.logsumexpfunkcję do zachowania precyzji. Jeśli masz jakąś funkcję lnyzwracającą logarytm, ymożesz to zrobić, np .:

z ssumy.misc import logsumexp
zaimportuj numpy jako np

xmin = 1e-15
xmax = 1e-5

# uzyskać wartości x rozmieszczone logarytmicznie
xvs = np.logspace (np.log10 (xmin), np.log10 (xmax), 10000)

# oceń swoją funkcję na xvs
lys = lny (xvs)

# wykonaj integrację reguły trapezu
deltas = np.log (np.diff (xvs))
logI = -np.log (2.) + logsumexp ([logsumexp (lys [: - 1] + delta), logsumexp (lys [1:] + delta)])

Wartość logIjest dziennikiem całki, którą chcesz.

Oczywiście to nie zadziała, jeśli musisz ustawić xmin = 0. Ale jeśli masz jakąś niezerową dodatnią dolną granicę całki, możesz po prostu bawić się liczbą punktów, xvsaby znaleźć liczbę, w której całka zbiega się.