Телескопы покупают здесь


A A A A Автор Тема: Численные методы решения задачи статической звезды.  (Прочитано 865 раз)

0 Пользователей и 1 Гость просматривают эту тему.

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
Численные методы решения задачи статической звезды.

Не знаю точно в какой раздел, поскольку особенно ничего дискуссионного тут пока нет,
и скорее тема будет интересна студентам , как пример для численного решения дифференциальных
уравнений. В качестве примера я взял систему дифференциальных линейных уравнений
1-го порядка для статической звезды из статьи Оппенгеймера-Волкова или что почти то же
Вайнберг , глава 11.

Точно не знаю , какой программой пользуются сейчас современные студенты и
научные сотрудники. Коды буду выкладывать согласно программе Maxima (wxMaxima) версия. 16.04.2
Это важно, не все версии работают корректно под 7-ку.
Неплохо было бы получить аналогичные коды в Mathlab, Mathematica..., если кто сможет это сделать.
Мой знакомый теоретик до сих пор работает в Фортране.
Заодно можно будет проверить те решения, которые я получил.
В основном использовал метод Рунге-Кутта 4 порядка, плагин уже встроен в эту версию.

Выяснилось, что численные решения в отличие от аналитических решений
имеют ряд особенностей и проблем. Среди них.

1. Задача обезрамеривания постоянных и параметров задачи
2. Проблема начальных данных или краевых условий и вопрос как их получить
3. Проблема нахождения особенности в решении
4. Проблема определения устойчивости решения численным способом.

Кроме них возникают (хотя у меня их невозникло) еще ряд проблем именно в численных методах.
- для сложных уравнений возможны осцилляции, расхождения решения...

Рассмотрю несколько примеров, поэтому тема не на 1 день.

1. Статическая звезда с однородной плотностью. Есть точное решение. Сравнение.
2. Звезда с ультрарелятивистким состоянием вещества ( есть точное решение).
3. Статическая звезда с разными политропами.
4. И если получится - Звезда с гибридным составом и Что такое Сверхмассивный Белый карлик?
5. Идеально было бы получить алгортм, куда , подставляя входные параметры, получил бы
составную звезду произвольной массы.

Итак, Оппенгеймер и Волков получили  основные уравнения для статической звезды (сборник
Альберт Эйнштейн и теория относительности).

\[ \frac{du}{dr}=4{\pi}{\epsilon}r^2 \quad(1) \]
\[ \frac{dp}{dr}=-\frac{p+\epsilon}{r(r-2u)}(4{\pi}r^3p+u) \quad(2) \]

Где \( \epsilon \) и \( p \) - плотность и давление , зависящее только от \( r \).

Здесь уже записаны уравнения , где безразмерны 2 константы : \( G=c=1  \)
 
ЗАДАЧА 1. Постоянная плотность.

\( \epsilon=const \)

Интегрируя уравнения (1) от \( 0 \) до \( r_b \) получим:
\[ u(r_b)-u(0)=r_g/2=(4/3){\pi}r_b^3\epsilon \]

Поскольку считают, что \( u(0)=0 \) , получим:

\[ \epsilon=\frac{3r_g}{8r_b^3{\pi}} \]

Смысл   функции \( u(r) \)  Это масса вещества плюс энергия гр. поля, заключенная  от \( r=0 \) до \( r \) . Имеется граничное условие , если звезда имеет границу \( r=r_b \):

\[ p(r=r_b)=0, u(r=r_b)=m=r_g/2 \]

\(  m \) - Шварцшильдовская масса. При данных константах размерность \( u \) - длина .
Таким образом у нас 2 произвольных параметра : масса звезды или  \(  r_g \)
и радиус  \(  r_b \)

Также у нас есть точное решение для давления , что позволит откалибровать программу.
Вайнберг, стр. 356, уравнение (11.6.4) А именно :
\[ p(r)=-\frac{\epsilon\left( \sqrt{-\frac{{{r}^{2}}\,{{r}_{g}}-{{{{r}_{b}}}^{3}}}{{{{{r}_{b}}}^{3}}}}-\sqrt{-\frac{{{r}_{g}}-{{r}_{b}}}{{{r}_{b}}}}\right) }{\sqrt{-\frac{{{r}^{2}}\,{{r}_{g}}-{{{{r}_{b}}}^{3}}}{{{{{r}_{b}}}^{3}}}}-3\sqrt{-\frac{{{r}_{g}}-{{r}_{b}}}{{{r}_{b}}}}}\quad(3) \]

Для примера возьму \( r_g=1 \) , \( r_b=10 \), шаг \( 1 \). от нуля до \( 10 \).
Особенность данных краевых условий в том, что у нас начинается расчеты не с начала \( r=0 \) ,  а с конца \( r=r_b \) В этом случае приходится Рунге Кутта задавать с отрицательным шагом \( -1 \).

res: rk([(4*r^2*%pi*epsilon),-(p+epsilon)*(u+4*%pi*r^3*p)/(r*(r-2*u))],[u,p],[r_g/2,0],[r,r_b,0,-r_b/20])

Коды Максима (приложил )

(кликните для показа/скрытия)



Как видно численное решение - красные точки - хорошо легли на синюю кривую точного решения. Что случилось не сразу , поэтому контроль важен.
« Последнее редактирование: 29 Июн 2018 [22:20:44] от ulitkanasklone »
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн Geen

  • *****
  • Сообщений: 12 210
  • Благодарностей: 200
  • Мне нравится этот форум!
    • Сообщения от Geen
Среди них.
Есть ещё одна проблема - точное решение может быть неустойчиво - для численных методов тогда мало шансов остаться влизи точного решения.
Если у тебя есть фонтан, заткни его, дай отдохнуть и фонтану.

А ещё мы любим обсуждать вкус устриц с теми кто их ел...

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
Есть ещё одна проблема - точное решение может быть неустойчиво - для численных методов тогда мало шансов остаться влизи точного решения.
Я думаю, мы столкнемся с такой проблемой. Уже во второй задаче.
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
Белый карлик я так понял получается с другим соотношением примерно \( r_b/r_g=6000 \).
Из этого численного решения по точкам видно, что непонятно  , где наступает сингулярность.
А из точного аналитического мы знаем это: \( r_b/r_g=9/8=1.125 \).
Это есть одна из проблем численных методов.

Зная точное решение, можно взять
\( r_b/r_g=1.13 \) и получить такой график:



По крайней мере рост давления в нуле виден.
На этом с задачей 1 можно завершить.
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
ЗАДАЧА 2. Ультрарелятивистский случай.

Ультрарелятивистское уравнение состояния
\[ \epsilon=3p \]
Тогда уравнения приобретают вид:

\[ \frac{du}{dr}=12{\pi}{p}r^2 \quad(4) \]
\[ \frac{dp}{dr}=-\frac{4p}{r(r-2u)}(4{\pi}r^3p+u) \quad(5) \]

Мне не удалось решить эту систему точно , но решение есть , оно приведено у Вайнберга (11.4.13) со ссылкой работу Мизнера с коллегами 60-х годов.

\[ p=\frac{1}{56{\pi}r^2}\quad(6) \]

Оно имеет особенность в \( r=0 \) и нет границы. И вот тут непонятно, как задавать
начальные данные и в какой точке. В принципе можно удалить особенность в нуле заменой
\[ x=pr^2 \]
тогда должно быть \( x=const=A \).
Должно быть: \( A=1/(56{\pi}) \)
Но изначально мы не знаем эту постоянную . Взяв точное решение (6) при \( r_0=1 \) , можно найти :

\( p_0=0.005684105110424833,  u_0=0.2142857142857143 \)

Тогда метод Рунге_Кутта дает хорошее соответствие.



Если постоянную \( A \) увеличить в 2 раза и взять начальные данные в той же точке \( r_0 \), то
получим такой график:



Если постоянную \( A \) уменьшить в 2 раза и взять начальные данные в той же точке \( r_0 \), то
получим такой график:



Видно , что расхождения есть  и это есть одна из проблем - задание начальных данных.

Коды Максима следующие (формулы уже подстроены под общий случай политропы):

(кликните для показа/скрытия)
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн Geen

  • *****
  • Сообщений: 12 210
  • Благодарностей: 200
  • Мне нравится этот форум!
    • Сообщения от Geen
В принципе можно удалить особенность в нуле заменой
x=pr2

тогда должно быть x=const=A.
Должно быть: A=1/(56π)
Но изначально мы не знаем эту постоянную . Взяв точное решение (6) при r0=1 , можно найти
Непонятно...
Если у тебя есть фонтан, заткни его, дай отдохнуть и фонтану.

А ещё мы любим обсуждать вкус устриц с теми кто их ел...

Оффлайн Geen

  • *****
  • Сообщений: 12 210
  • Благодарностей: 200
  • Мне нравится этот форум!
    • Сообщения от Geen
И вот тут непонятно, как задавать
начальные данные и в какой точке.
Их нет.
Вообще, фотонный газ, что в этом примере, что в случае всей вселенной, не даёт возможность выбрать реперы/масштабы...
Если у тебя есть фонтан, заткни его, дай отдохнуть и фонтану.

А ещё мы любим обсуждать вкус устриц с теми кто их ел...

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
Непонятно...
Почитайте стр. 344 у Вайнберга. Мне тоже его рассуждения непонятны. В данном случае я имел в виду, что функция \( r^2p \) по крайней мере интегрируется без особенности. По смыслу \( r^2p \) это сила действующая на сферу.
И полная масса \( u(r) \) интегрируется.
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн Geen

  • *****
  • Сообщений: 12 210
  • Благодарностей: 200
  • Мне нравится этот форум!
    • Сообщения от Geen
Мне тоже его рассуждения непонятны.
Но не понятно что именно Вы пишете (и делаете) начиная с процитированного фрагмента...
Если у тебя есть фонтан, заткни его, дай отдохнуть и фонтану.

А ещё мы любим обсуждать вкус устриц с теми кто их ел...

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
Но не понятно что именно Вы пишете (и делаете) начиная с процитированного фрагмента...
Мне в таком (изначальном виде)  не удалось решить систему точно.
Такой заменой \( x=r^2p \) второе выражение несколько упрощается. Возможно можно догадаться о точном решении. То есть догадаться, что \( x=const \) . И в новых функциях можно пытаться численно решать систему начиная с \( r=0 \) . Я правда этого не сделал.
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
3. ЗАДАЧА 3. Звезда из вещества с состоянием политропы.

Тут меня ждала неудача.
Можно взять результаты отсюда:
http://www.astronet.ru/db/msg/1252779/42.html
Для нерелятивистского случая , когда \( r_b>>r_g \)
система уравнений упрощается:
\[ \frac{dp}{dr}=-\frac{{\epsilon}u(r)}{2r^2} \]
\[ \frac{du}{dr}=4{\pi}r^2{\epsilon} \]
\[ p=K_{\gamma}{\epsilon}^{\gamma} \]

Или последнее можно записать \( \epsilon=K_1p^s \), \( s=1/\gamma \)

Взял как пример -  Белый карлик с массой Солнца и с размером \( r_b/r_g=5000  \),
политропу \( \gamma=5/3 \), \( K_1=10^{-8} \).
Граничные условия я ставил как в первой задаче : на границе давление ноль \( p(r=r_b)=0 \) и \( u(r=r_b)=r_g/2 \) ,
численный метод Рунге_Кутта не дал вменяемый результат, а только тривиальный \( p=0 \) и \( u=r_g/2 \).

g:4*%pi*r^2*epsilon;
f:((-p-p^(s))*(u+4*%pi*p*r^3))/(r*(r-2*u))
res: rk([g,f],[u,p],[r_g/2,0],[r,r_b,0,-r_b/n])

Либо я ставлю нереальный граничные условия, либо Рунге-Кутт не работает
в данном случае. Пока не могу разобраться. Если ничего
не придумаю или никто не подскажет, то перейду к 4 примеру, а потом вернусь к этому.


 
« Последнее редактирование: 03 Июл 2018 [12:24:37] от ulitkanasklone »
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн VladTK

  • *****
  • Сообщений: 2 182
  • Благодарностей: 61
  • Через тернии к звездам
    • Сообщения от VladTK
Возможно Вам будет интересно https://arxiv.org/abs/1806.05719
Celestron C6-N

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
Возможно Вам будет интересно https://arxiv.org/abs/1806.05719
Спасибо попробую разобраться. Но я так вижу сходу, они просто взяли \( \gamma(r) \) переменную. Она должна быть от 4/3 до 5/3 . Но как именно она зависит непонятно.
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html

Оффлайн ulitkanaskloneАвтор темы

  • *****
  • Сообщений: 8 645
  • Благодарностей: 280
  • тихо,тихо ползи улитка по склону Фудзи..
    • Сообщения от ulitkanasklone
    • СПРАВОЧНАЯ ИНФОРМАЦИЯ ПО ОТО, СТАТЬИ И НЕКОТОРЫЕ ЗАДАЧИ
Возможно Вам будет интересно
Все таки тут дело в математическом методе.

Нашел хорошую методичку Иванова по звездам и политропам и взял его дифференциальный уравнения на стр. 157.
http://lnfm1.sai.msu.ru/~rastor/Study/Ivanov_physicsofstars.pdf
Взял готовые уравнения со стр. 157:
\[ \frac{dp}{dx}=-q\sigma/x^2 , \quad \frac{dq}{dx}=x^2\sigma , \quad p=\frac{c}{{(4{\pi)}}^{1/n}}{\sigma}^{1+1/n} \]

Где все обезразмерено и для \( \gamma=5/3  \)  соответствует \( n=3/2 \) и \( c=0.42 \).

Краевые условия  : \(  p'=0 (x=0) \) и \( q=1 , p=0 (x=1) \)

Рунге Кутт дает тот же невразумительный результат, все нули. Значит надо решать как-то все по другому , не с использованием данных Коши в конечно точке. Либо сводить все к уравнению второго порядка, хотя мне  это не советовали.

Точнее , если на границе ставить давление не ноль , то можно получить профиль для давления и массы (плотности). Это бы означало, что обязательно у звезды есть атмоосфера и нет четкой границы.
« Последнее редактирование: 09 Июл 2018 [13:54:12] от ulitkanasklone »
Мне нравится этот форум
моя страница:
http://антониум.рф/ОТО/GT.html