Использование простого метода Эйлера для решения EoS звезд (Mathematica) - PullRequest
Новая
       27

Использование простого метода Эйлера для решения EoS звезд (Mathematica)

2 голосов
/ 17 апреля 2020

Я попытался перевести код Matlab: 

clear
clc;

%Use a resolution of dr = 10^-N
N = 4;

%G = GAMMA, adiabatic parameter
G = 4/3; % high density/extreme relativistic
%G = 5/3; % low density/nonrelativistic

r = 0;
dr = 10^-N;
m(1) = 0;
r(1) = 0;
rho(1) = 0.1;
i = 1;
    while (i>0)
        i=i+1;
        r(i) = r(i-1)+dr;
        Gamma(i) = rho(i-1)^(G-1) / (3*sqrt(2)) ;%(rho(i-1))^(2/3) / (3*(1+rho(i-1)^(2/3) )^(1/2) );

        drho = (m(i-1)*rho(i-1)) / (Gamma(i) * r(i)^2)*dr;
        if (drho > rho(i-1))
            break;
        else
            rho(i) = rho(i-1) - drho;
        end
        m(i) = m(i-1) + (r(i)^2 * rho(i)) * dr;
    end

    Mass = m(i-1)
    Rad = r(i-1)

plot(r(1:i-1),m(1:i-1))
xlabel('Log10(r/[m])')
ylabel('M(r) (M_{sol})')

В Mathematica, что у меня есть: 

EoS[n0_] :=
 Module[{n = n0},
  N1 = 4;
  G1 = 4/3;
  G2 = 5/3;
  rho[n_] := 0.1 + n;
  m1[n_] := n;
  dr = 10^(-N1);
  r[n_] := n;

  Gamma1[n_] := ((rho[n - 1])^(G1 - 1))/(3*Sqrt[2]);
  drho = ((m1[n - 1]*rho[n - 1]))/(Gamma1[n - 1]*r[n - 1]^2)*dr;
  n = 1;
  While[n > 0, Gamma1[n] && drho; n++]
    If[drho > rho[-1], Break[], rho[n] = rho[-1] - drho]
    m1[n] = m1[-1] + ((r[n]^2)*rho[n])*dr;

  Mass = m1[n - 1];
  Rad = r[n - 1];
  ParametricPlot[{m[n], r[n]}, {n, 1, n - 1}, 
   AxesLabel -> {"Log10(r/[m])", "M(r)(M_{sol})"}]
  ]

Моя проблема в том, что я получаю сложные решения, когда я, очевидно, не должен «т. Я абсолютный новичок в математике, особенно с модулями и самим методом Эйлера. Любая обратная связь / помощь будет высоко ценится. Я предполагаю, что мои приращения неверны, я действительно не знаю.

...