Całka w przestrzeni log-log

10

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źć stfa(r)

fa(r)=0ry(x)rex

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

DilithiumMatrix
źródło
Znalazłem ten link ( my.originlab.com/forum/topic.asp?TOPIC_ID=1251 ), który wydaje się podążać tą samą drogą, którą normalnie bym poszedł: obliczanie nachylenia i przechwytywanie w przestrzeni dziennika. Następnie przekonwertuj na przestrzeń lin-lin, zintegruj i oceń.
MrMas

Odpowiedzi:

6

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

F(r)=log(r)exp(a+b)da

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

Damaszek Stal
źródło
Dzięki za twoją odpowiedź! Ale myślę, że to skutecznie wykonuje całkę w przestrzeni liniowej. Być może jednak proszę o coś niemożliwego ...
DilithiumMatrix
2
Nie, to robi całkę w przestrzeni dziennika. Dyskretyzującrezajest równej wielkości w przestrzeni dziennika, a nie przestrzeni liniowej.
Damascus Steel
1
@DilithiumMatrix ma rację: dyskretyzacja x wartości są w przestrzeni dziennika, ale interpolacja ywartości zachodzą w przestrzeni liniowej. Tak więc, jeśli użyjesz reguły trapezoidalnej, funkcja, która jest skutecznie zintegrowana, jest fragmentarycznie liniowa na wykresie z logarytmiczną osią x i liniową osią y.
burnpanck 30.04.2018
3

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

F(r)=0ry(x)dx
myślisz coś takiego
F=integrate(y,x)
gdzie F=[F1,...,Fn] jest wektorem próbkującym całkę na siatce x.

Jednak nie masz próbek x i y, ale raczej masz próbki x^=log(x) i y^=log(y).

Oczywiście najprostsze podejście byłoby

F=integrate(exp(y^),exp(x^)),
ale byłoby to podatne na błędy, ponieważ y(x) mimo to nie jest gładka y^(x^) jest.

Teraz reguła trapezoidalna zasadniczo zakłada twój wkłady(x)jest fragmentarycznie liniowy. Zatem proste uogólnienie byłoby takie założeniey^(x^) jest fragmentarycznie liniowy.

W tym przypadku definiowanie ΔFk=Fk+1Fk, ty masz

ΔFk=xkxk+1y(x)dx=x^kx^k+1ey^(x^)ex^dx^=x^kx^k+1y~(x^)dx^

Następnie definiowanie t=(x^x^k)/Δx^k, ty masz

y^k+ty^k+tΔy^k
i y~(t)aebt, z a=ey^k+x^k i b=Δy^k+Δx^k.

Tak więc całka staje się

ΔFkaΔx^01ebtdt=aΔx^eb1b

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 y^(x^) dyskretnie dyskretnie x^ siatka, z F(x^1)=0.)

GeoMatt22
źródło
Dziękuję za twoją (bardzo jasną) odpowiedź, ale jak już powiedziałem w odpowiedzi na @DamascusSteel --- Myślę, że to po prostu odwraca całkę do liniowo-liniowej przestrzeni i traci zalety przestrzeni logarytmicznej.
DilithiumMatrix
1
@DilithiumMatrix: To nie jest ta sama odpowiedź, co DamascusSteel. Zauważ, że zastosowanie reguły trapezoidalnej do odpowiedzi DamascusSteel nie dałobyexp(b)1bczynnik.
burnpanck
3

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(xi,yi) i (xi+1,yi+1) przy założeniu, że y=Cixni w przedziale czasowym i, otrzymujesz ni=ln(yi/yi+1)/ln(xi/xi+1) and Ci=ln(yi)niln(xi). The contribution to the integral from interval i is then

ΔFi=xixi+1Cixnidx={Cini+1(xi+1ni+1xini+1),ni1Ci(lnxi+1lnxi),ni=1,
where you obviously need some tolerance for identifying the special case ni=1 in your implementation.
Stefan B. Lindström
źródło
3

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: wprowadź opis zdjęcia tutaj

log(y1)=m1log(x1)+n1
log(y2)=m1log(x1)+n1

where m1 is the slope of the line and n1 is its y-intercept.

By subtracting the two, one can find:

m1=log(y1)log(y2)log(x1)log(x2)

i od podstawienia:

n1=log(y1)-m1log(x1)

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)xmmin

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

x1x2)y(x)rex=min1m1+1(x2)m1+1-x1m1+1),dla m-1

i

x1x2)y(x)rex=min1logx2)x1,dla 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.

Elena Pascal
źródło
This is wonderful @elenapascal, this has been bothering me for over 3 years now, and I think this is (or is very close to) the solution. I don't quite following your last relation, I don't think the integral over y equal to the log(x2/x1)
DilithiumMatrix
In particular, if I take the log of the integral on the left-hand-side, then I get a similar term to the right hand side, but with log([x_2/x_1]^{m_1+1} + 1), i.e. there is an additional +1 in the argument of the log
DilithiumMatrix
Bardzo mnie to dzisiaj niepokoiło, dopiero po napisaniu tego zdałem sobie sprawę, że @Stefan opublikował tę samą odpowiedź. Dla m = -1 wystarczy zastąpić to w definicji y: y (x) = e ^ n / x. To daje dzienniki. Nie jestem pewien, czy śledzę twój drugi post
Elena Pascal
Właśnie zdałem sobie sprawę z tego samego, ale nie do końca zrozumiałem, dopóki nie przeczytałem twojego wyjaśnienia
DilithiumMatrix
1

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ę.

Matt Pitkin
źródło