Численное значение интеграла

mtk79

Дано некоторое малое число h, например, h=10^{-4}
далее вводятся следущие величины как функции угла t
[math] $$ \xi=1-\sqrt{1-h^2}\cos t;\,    B=\frac{h \sin t}{\xi}    $$ [/math]
мне нужно численное значение интеграла (или главную степень лорановского разложения по h)
[math]$$ \int_0^{\pi}  \frac{\sin^3 t}{(\xi+B)^2} dt    $$ [/math]
например, при h=10^{-4}.
Я это пытался провести в Maple 7 и Maple 12 — но они либо виснут, либо выдают эпические фейлы. Если кто может провести у себя или в других матпакетах — то Ваше имя я готов буду прославить на доступных мне просторах форума.локал (т.е. там, где для меня еще не выписали рестрикт).
а сам я уже отвык от вычисления угловых интегралов (стыдно за МГУ)
Чтобы не переписывать, код можно увидеть и скопировать, отключив картинку.

Andrey56

Mathematica говорит 1.98778

mtk79

спасибо. А задачу нахождения главного (как Вы показали, нулевого) члена разложения по степеням h может? Т.е.можно ли извлечь аналитическое число, которое публиковать гораздо приятнее?
Исходя из Вашего рез-та, интуиция подсказывает, что это 2+с*sqrt{h} или что-то в этом духе

Andrey56

Пардон, ответ 9.88529
Квадрат забыл поставить.

Andrey56

Это думать надо, а я не умею )

mtk79

Все равно спасибо. Правда, моя квантовая интуиция потерпела мощный фейл после уточнения. Под каким именем Вас прославлять в доступных мне чатах (как я понимаю, прочтение Вашего ника по-гишпански как "Иипок" не катит)?

lenmas

Судя по всему, если разложение в ряд Лорана и есть, то имеет полюс не более первого порядка (в знаменателе сумма взаимно обратных величин оценивается через
[math]  $  \sqrt{h\sin t}  $  [/math] снизу) :)
Да, кстати, числовые значения интегралов вычиляются в Matlab'е.
И еще одно кстати, интеграл должен не так уж плохо браться, если под дифференциал занести один из синусов.
P.S. Не, по поводу последнего я че-то загнался. Синус-то один останется :grin:

seregaohota

Напиши в Maple 7 кроме первой строчки (> interface(prettyprint=0) чтобы видеть нормально.
При h=1e-4 действительно выдаёт то значение, которое и математика считает, при других h естественно другой ответ.
> interface(prettyprint=0);
> assume(h>0);
> assume(epsilon>0);
> xi:=1-sqrt(1-h^2)*cos(t); B:= h*sin(t)/xi;
xi := 1-(1-h^2)^(1/2)*cos(t)
B := h*sin(t)/(1-(1-h^2)^(1/2)*cos(t
> J:=Int(sin(t)^3/(xi+B)^2,t=epsilon..Pi);
J := Int(sin(t)^3/(1-(1-h^2)^(1/2)*cos(t)+h*sin(t)/(1-(1-h^2)^(1/2)*cos(t^2,t = epsilon .. Pi)
> S:=series(J,h,2);
S := series-2*(2*ln(tan(1/2*epsilon+2*ln(tan(1/2*epsilon*tan(1/2*epsilon)^2-ln(1+tan(1/2*epsilon)^2)-ln(1+tan(1/2*epsilon)^2)*tan(1/2*epsilon)^2+1)/(1+tan(1/2*epsilon)^2+(-2/3*(3*Pi*tan(1/2*epsilon)^3+2-6*tan(1/2*epsilon)^2-6*arctan(tan(1/2*epsilon*tan(1/2*epsilon)^3)/tan(1/2*epsilon)^3)*h+O(h^2h,2)
> evalf(subs(epsilon=1e-4,S;
> evalf(subs({epsilon=1e-4,h=1e-4},J;
series(37.61395018-.1066666656e14*h+O(h^2h,2)
9.885286008
> evalf(subs(epsilon=1e-4,S;
> evalf(subs({epsilon=1e-4,h=1e-5},J;
series(37.61395018-.1066666656e14*h+O(h^2h,2)
12.94451076
> S0:=coeff(S,h,0);
S0 := -2*(2*ln(tan(1/2*epsilon+2*ln(tan(1/2*epsilon*tan(1/2*epsilon)^2-ln(1+tan(1/2*epsilon)^2)-ln(1+tan(1/2*epsilon)^2)*tan(1/2*epsilon)^2+1)/(1+tan(1/2*epsilon)^2)
> limit(S0,epsilon=0,right);
infinity
>

lenmas

Ты лучше попробуй посчитать предел после умножения на h или на корень из h :)
Подобрать степень h, по которой он уходит на бесконечность.
В wxMaxima у меня тоже такие же ответы получились :)

lenmas

Чего-то я прикинул по минус степеням десятки первые четыре или пять значений от -4-ой степени до -8-ой. В общем расстояния между ними стабилизируются, то-есть очень похоже на c*log(1/h причем коэффициент при логарифме у меня на каждом шаге все ближе и ближе к 4/3 становится. Не знаю, насколько это правда, но автор треда вообще-то мог бы уж как-нибудь и взять интеграл, интеграл-то берущийся. Только побольше бумаги запасти надо, и большой по размеру :grin: И использовать вместо sqrt(1-h^2) и h какие-нибудь нормальные обозначения, ну cos(h) и sin(h) что ли. Неужели для статьи не пойти на такую жертву? ;)

mtk79

Если бы весь интеграл был I(h) \sim 4/3*log(1/h) +O(1) — то как при 1/h=10000 получается 1.88?

seregaohota

У тебя же вроде Maple 7 доступна ты писал. Не повторил в Maple 7 что ли, что я с выводом на экран синенькими буковками в текстовом виде тут в запостил так как мне влом было картинки в форуме набирать и упрощать кое-что "врукопашную"? :) Если ты повторишь в Maple 7 всё что в том моём посте было или загонишь в LaTeX тот мой пост, то увидишь в нормальном математическом виде.
Я ввёл второй малый параметр [math]$\varepsilon \to +0$[/math], интеграл не от 0, а от эпсилон считал.
Короче на пальцах что выдает Maple 7.
Я посмотрел асимптотику - 2 члена разложения в ряд интеграла J по степеням h приведены, и его предел при h->0 посчитан [math]$S0 \to \infty$[/math], асимптотика S0 по [math]$\varepsilon \to +0$[/math]:
[math]$S0 \sim -4\ln\left( \tg \frac{\varepsilon}{2} \right)$[/math]. Дело в том я так понимаю, что при фиксированном h получается конкретный численный предел по [math]$\varepsilon \to +0$[/math], так как [math]$J = S0 + h S1 + O(h^2)$[/math] и в S0 и в S1 особенности по эпсилон взаимно компенсируются. При h=0 этого не происходит, т.к. S0 уходит в бесконечность.
PS Если Maple 7 не врёт, иногда бывает...
PPS При epsilon=1e-4 численные значения и ряда, и интеграла я там привожу при h=1e-4 и 1e-5, это 37.61395018-.1066666656e14*h+O(h^2) = 9.885286008 и 12.94451076 соответственно.
Если взять epsilon=1e-7 например, то при h=1e-4 всё так же, а при h=1e-5
65.24497132 - .1066666667e23*h + O(h^2) = 12.94451076
И S0 и S1 уходят в бесконечность, но не по h, а по эпсилон, и при фиксированном h взаимно компенсируются, численные значения в ппределе от эпсилон при фиксированном h не меняются.

lenmas

Там же при h=10^{-4} получается 9.88... :confused: Ну и никто не отменял достаточно маленького слагаемого в O(1). Я смотрел при достаточно маленьких h, 10^{-4}, 10^{-5}, 10^{-6}, 10^{-7} и 10^{-8}. Так вот расстояния между значениями интеграла при них почти постоянны. А константу в O(1) я не прикидывал, хотя можно было бы тоже по значениям интегралов, но я уже закрыл программу давно. В общем, сейчас посмотрел, константа стремится к -2.409 примерно, то-есть асимптотика типа 4/3 log(1/h)-2.409+O(h). Но я бы все-таки честно вычислил интеграл, и посмотрел аналитическое выражение для интеграла. Может, там просто число очень похожее на 4/3, а реально стоят какие-нибудь корни. :grin:
P.S. Да, если ты не понял, все логарифмы имеются в виду натуральные, конечно.

mtk79

Да, спасибо. Пост я видел, а вот воспроизвести не мог, ибо для Мапла-7 нужна была перезагрузка компа с другого диска, я ждал, когда все юзеры освободят комп.

mtk79

Я что-то тупанул. Успокоился, что 1.98 и уточненное уже прочитал как 1.89, а не 9.

lenmas

Сегодня поигрался с Maple 11. В общем, гипотеза про 4/3 подтверждается. А вот по поводу константы в O(1) что-то не уверен. Не мог бы кто-нибудь с мощным компом зааплодить сюда вывод мэпла на запрос
[math]  $$  int\Bigl(\frac{\sin(x)^3}{(1-(1-b^2/2)\cos(x)+b\sin(x)/(1-(1-b^2/2)\cos(x^2},x=0..\mathrm{Pi}\Bigr);  $$  [/math]
?
А то уходил с работы, так и не дождался, когда он выплюнет ответ :)

lenmas

Отпишусь-ка последний раз в этом треде. :grin:
Судя по маплу, асимптотика имеет вид
[math]  $$  \frac43\ln\Bigl(\frac2h\Bigr)-\frac{10}3+O(h\quad h\to0.  $$  [/math]
Оставить комментарий
Имя или ник:
Комментарий: