Численные методы решения задачи статической звезды.Не знаю точно в какой раздел, поскольку особенно ничего дискуссионного тут пока нет,
и скорее тема будет интересна студентам , как пример для численного решения дифференциальных
уравнений. В качестве примера я взял систему дифференциальных линейных уравнений
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])
Коды Максима (приложил )
/* [wxMaxima: input start ] */
r_g:1; r_b:10;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
float(9/8);
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
epsilon:3*r_g/(8*%pi*r_b^3);
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
'diff(u,r)=4*%pi*r^2*epsilon;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
'diff(p,r)=-(p+epsilon)*(u+4*%pi*r^3*p)/(r*(r-2*u));
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
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])$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
res;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
'p1=factor((epsilon*sqrt(1-(r^2*r_g)/r_b^3)-epsilon*sqrt(1-r_g/r_b))/(3*sqrt(1-r_g/r_b)-sqrt(1-(r^2*r_g)/r_b^3)));
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
p1:factor((epsilon*sqrt(1-(r^2*r_g)/r_b^3)-epsilon*sqrt(1-r_g/r_b))/(3*sqrt(1-r_g/r_b)-sqrt(1-(r^2*r_g)/r_b^3)));
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
solp:makelist([y[1],y[3]],y,res);
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
solu:makelist([y[1],y[2]],y,res);
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
plot2d ([[discrete, solp],p1],[r,0.01,r_b],
[style,points,lines], [point_type,diamond], [color,red,blue],
[legend, "численно", "точное"],
[xlabel, "r/r_g"],
[ylabel, "давление p"])$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
plot2d ([[discrete, solu],
[style,points], [point_type,diamond], [color,red],
[xlabel, "r/r_g"],
[ylabel, "давление p"])$
/* [wxMaxima: input end ] */
Как видно численное решение - красные точки - хорошо легли на синюю кривую точного решения. Что случилось не сразу , поэтому контроль важен.