Интеграл в пространстве журнала

10

Я работаю с функциями, которые, в общем, гораздо более плавные и лучше ведут себя в пространстве журналов регистрации - так что именно здесь я выполняю интерполяцию / экстраполяцию и т. Д., И это работает очень хорошо. Есть ли способ интегрировать эти числовые функции в пространстве журнала?

т.е. я надеюсь использовать какое-то простое трапециевидное правило для выполнения кумулятивного интеграла (например, в python, use scipy.integrate.cumtrapz), чтобы найти некоторый stF(р)

F(р)знак равно0рY(Икс)dИкс

Но я надеюсь использовать значения и вместо и (когда это возможно).Lог(Y)Lог(Икс)YИкс

DilithiumMatrix
источник
Я нашел эту ссылку ( my.originlab.com/forum/topic.asp?TOPIC_ID=1251 ), которая, кажется, движется так же, как обычно: вычислять наклон и перехват в пространстве журнала-журнала. Затем преобразуйте в пространство lin-lin, интегрируйте и оцените.
MrMas

Ответы:

6

Вы можете просто изменить переменные. Установка , . Интеграл становитсяa=log(x)b(a)=log(y(x))

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

Вы должны быть немного осторожнее, потому что вы интегрируете из . Что именно вы должны сделать, будет зависеть от того, как выглядит .y(x)

Дамасская сталь
источник
Спасибо за ваш ответ! Но я думаю, что это просто выполнение интеграла в линейном пространстве. Возможно, я прошу что-то невозможное, однако ...
DilithiumMatrix
2
Нет, это делает интеграл в лог-пространстве. При дискретизации одинаково измеряется в лог-пространстве, а не в линейном пространстве. da
Дамаск Сталь
1
@DilithiumMatrix прав: дискретизация значений находится в лог-пространстве, но интерполяция значений происходит в линейном пространстве. Таким образом, если использовать правило трапеции, эффективно интегрированная функция является кусочно-линейной на графике с логарифмической осью X и линейной осью Y. xY
Burnpanck
3

Я не использую python, но если я правильно понимаю, то с помощью вы думаете о чем-то вроде где - вектор, в котором выполняется выборка интеграла по сетке .

F(r)=0ry(x)dx
F=integrate(y,x)
F=[F1,...,Fn]x

Однако у вас нет образцов и , а есть образцы и .xyx^=log(x)y^=log(y)

Конечно, самый простой подход будет

F=inTегрaTе(ехр(Y^),ехр(Икс^)),
но это было бы подвержено ошибкам, потому что Y(Икс) не гладко, хотя Y^(Икс^) является.

Теперь трапециевидное правило по существу предполагает ваш вкладy(x)кусочно-линейный Таким образом, простое обобщение будет для вас предположить, чтоy^(x^) кусочно-линейный

В этом случае, определяя ΔFk=Fk+1Fk, у тебя есть

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

Затем, определяя t=(x^x^k)/Δx^k, у тебя есть

y^k+ty^k+tΔy^k
а также y~(t)aebt, с a=ey^k+x^k а также b=Δy^k+Δx^k,

Таким образом, интеграл становится

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

В Matlab это будет выглядеть примерно так

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]);

Надеюсь это поможет!

(Редактировать: мой ответ по сути такой же, как и гораздо более краткий ответ, который Дамаск Сталь дал, когда я печатал. Единственное отличие состоит в том, что я пытался дать конкретное решение для случая, когда «конкретный y(x)"кусочно-линейный y^(x^) дискретизирован по дискретному x^ сетка, с F(x^1)=0.)

GeoMatt22
источник
Спасибо за ваш (очень четкий) ответ, но, как я только что сказал в ответ на @DamascusSteel - я думаю, что это просто возврат интеграла к линейно-линейному пространству и потеря преимуществ лог-пространства.
DilithiumMatrix
1
@DilithiumMatrix: это не тот ответ, который дал DamascusSteel. Обратите внимание, что применение правила трапеции к ответу DamascusSteel не дастехр(б)-1бфактор.
Burnpanck
3

Если функция выглядит гладко на графике log-log, вы можете интерполировать, используя степенной закон на каждом интервале (степенные законы, конечно, линейны в log-log). Таким образом, между точками(xi,yi) а также (xi+1,yi+1) при условии, что y=Cixni в пределах интервала i, вы получаете ni=ln(yi/yi+1)/ln(xi/xi+1) а также Ci=ln(yi)niln(xi), Вклад в интеграл от интервалаi затем

ΔFi=xixi+1Cixnidx={Cini+1(xi+1ni+1xini+1),ni1Ci(lnxi+1lnxi),ni=1,
где вам, очевидно, нужна некоторая терпимость для выявления особого случая ni=1 в вашей реализации.
Стефан Б. Линдстрем
источник
3

Я думаю, что есть некоторая путаница с изменением переменных в некоторых из предыдущих ответов, а также с некоторыми ошибками. Интеграл лог-функции - это не лог интеграла. Я думаю, что в целом сложно выписать интеграл функции, зная интеграл ее логарифма. Если кто-нибудь знает, как это сделать, мне было бы интересно.

Между тем, решение @ Stefan, представленное выше, является способом обойти интеграцию функции в пространстве log-log. Отправной точкой является то, что функция, с которой вы работаете, является линейной в пространстве log-log для достаточно маленьких сегментов.

Затем можно написать уравнение линии в конечных точках сегмента: введите описание изображения здесь

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

где m1 это наклон линии и n1 это у-перехват.

Вычитая два, можно найти:

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

and from substitution:

n1=log(y1)m1log(x1)

If in the log-log space the equation of a segment is close to a line then in normal (linear) space the equation of the segment is close to an exponential:

y(x)xmen

If we have a analytical formulation for this segment it is easy to integrate:

x1x2y(x)dx=en1m1+1(x2m1+1x1m1+1),for m1

and

x1x2y(x)dx=en1logx2x1,for m=1

This feels a bit like cheating, but this is sampling in log-log space such that we can approximate the function in the linear space to an exponential with parameters derived from the log-log space.

Elena Pascal
источник
Это замечательно @elenapascal, это беспокоит меня уже более 3 лет, и я думаю, что это (или очень близко) решение. Я не совсем после вашего последнего отношения, я не думаю , что интеграл по у равен к бревну (x2 / x1)
DilithiumMatrix
В частности, если я возьму лог интеграла с левой стороны, то получу аналогичный термин с правой стороны, но с log([x_2/x_1]^{m_1+1} + 1), т. Е. В аргументе журнала есть дополнительный +1
DilithiumMatrix
Меня это очень беспокоило и сегодня, только после того, как я написал это, я понял, что @Stefan опубликовал тот же ответ. Для m = -1 вы просто замените это в определении y: y (x) = e ^ n / x. Это дает журналы. Я не уверен, что следую за вашим вторым постом
Елена Паскаль
Я просто понял , то же самое, но я не до конца понял , пока я не прочитал ваше объяснение
DilithiumMatrix
1

Решение, которое я использую, является в основном реализацией правила трапеции и использует scipy.misc.logsumexpфункцию для поддержания точности. Если у вас есть какая-то функция, lnyкоторая возвращает логарифм, yвы можете сделать, например:

из scipy.misc импорт logsumexp
импортировать numpy как np

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

# получить значения х логарифмически
xvs = np.logspace (np.log10 (xmin), np.log10 (xmax), 10000)

# оцените свою функцию на xvs
lys = lny (xvs)

# выполнить интеграцию правила трапеции
deltas = np.log (np.diff (xvs))
logI = -np.log (2.) + logsumexp ([logsumexp (lys [: - 1] + deltas), logsumexp (lys [1:] + deltas)])

Значение logI- это журнал интеграла, который вы хотите.

Очевидно, это не будет работать, если вам нужно установить xmin = 0. Но если у вас есть некоторый ненулевой положительный нижний предел для интеграла, вы можете просто поиграть с количеством точек, xvsчтобы найти число, где интеграл сходится.

Мэтт Питкин
источник