1
$\begingroup$

My naive attempt:

ClearAll["Global`*"]
Md = 10^(-9); 
P = 10; 
R = 10^4; 
α = 10^(-2); 
ϵ = 10^(-4); 
γ = 10^(-2); 
ke = 0.02*(1 + 0.6625); 
k0 = 5*10^20; 
σ = 5.67*10^-8; 
Rg = 8315; 
c = 3*10^8; 
G = 6.67/10^11; 
M = 2.8*10^30; 

Ωk[r_] := Sqrt[(G*M)/r^3]; 
μ = Md/(3*Pi); 
κ = (27*ke)/(2*σ) Rg/μ; 
Co[r_] := 1; 
β[r_] := 0; 

Do[Σ[
   r_] := κ^(-1/5) α^(-4/5) μ^(3/
      5) Ωk[r]^(2/5)*Co[r]^(-1/5);
 h[r_] := (κ α Σ[r]^2 Ωk[
       r]^(-5) Co[r])^(1/6);
 T[r_] = 1/2 μ/Rg (Ωk[r]^2 h[r]^2)/(1 + β[r]);
 Kkr[r_] := k0 Σ[r]/h[r]*T[r]^(-7/2);
 β[r_] := μ/Rg (4 σ)/(3 c)
    T[r]^3/Σ[r] h[r];

 (*Co[r_] := (1 + β[r])^4*(1 + Kkr[r]/ke)*), {2}]

     Plot[Σ[r],{r,10^4, 10^10}]

     Plot[h[r],{r,10^4, 10^10}]

     Plot[T[r],{r,10^4, 10^10}]

to solve iteratively the system of equations:

$$\begin{array} &\mu=\dfrac{M_d}{3\pi}\\ &\rho(r) = \Sigma(r)\, h(r)\\ &\Omega_K(r)= \sqrt{\dfrac{G M}{r^3}}\\ &\beta(r) = \dfrac{\mu}{Rg}\, \dfrac{4\sigma}{3c}\, \dfrac{T^3(r)}{\rho(r)}\\ &T(r) =\dfrac{1}{2} \, \Omega_k^2(r)\, h^2(r)\, \dfrac{1}{1+\beta(r)}\, \dfrac{\mu}{Rg} \\ &h^6(r)=\kappa \alpha \Sigma^2(r) \, \Omega_k(r)^{-5}C(r)\\ &\Sigma(r) = k^{-1/5}\,\alpha^{-4/5}\mu^{3/5}\, \Omega_k^{2/5}(r)\, C^{-1/5}(r)\\ &C(r)=(1-\beta(r))^4 \left( 1+\dfrac{Kkr}{ke} \right)\\ &Kkr(r)=ko \rho(r) \, T^{-7/2}\\ \end{array}$$

fails. How can I solve it?

The iteration procedure goes like this:

  • Impose $C=1$ and $\beta(r)=0$
  • Solve for $\Sigma$
  • Solve for $T$ and $h$
  • Calculate $C$ and $\beta$ and $Kkr(r)$
  • Solve for $\Sigma$
  • ...
$\endgroup$
  • $\begingroup$ Does FindRoot not find the solution when you give it decent starting values? $\endgroup$ – Sjoerd Smit Jan 31 at 14:31
  • $\begingroup$ @SjoerdSmit I need to plot the functions $\Sigma$ $h$ and $T$. $\endgroup$ – mattiav27 Jan 31 at 14:40
  • $\begingroup$ It should possible to use FindRoot for a range of values of r, solving the equations for each individual value, right? In fact, once you have a solution for 1 value of r, you could use that solution as an initial guess for neighbouring values of r. $\endgroup$ – Sjoerd Smit Jan 31 at 15:09
  • $\begingroup$ Code and equations do not match. $\endgroup$ – Alex Trounev Jan 31 at 15:09
  • $\begingroup$ @AlexTrounev I ìt should be fine now $\endgroup$ – mattiav27 Jan 31 at 15:13
3
$\begingroup$

We can use the method of the false transient:

ClearAll["Global`*"]
Md = 10^(-9);
P = 10;
R = 10^4;
α = 10^(-2);
ϵ = 10^(-4);
γ = 10^(-2);
ke = 0.02*(1 + 0.6625);
k0 = 5*10^20;
σ = 5.67*10^-8;
Rg = 8315;
c = 3*10^8;
G = 6.67/10^11;
M = 2.8*10^30;

Ωk[r_] := Sqrt[(G*M)/r^3];
μ = Md/(3*Pi);
κ = (27*ke)/(2*σ) Rg/μ;
ic = {Co[0] == 1, β[0] == 
    0, Σ[
     0] == (μ^(3/5)*Ωk[
         r]^(2/5)) (κ^(-1/5)*α^(-4/5)*Co[0]^(-1/5)),
   h[0] == (κ*α*Σ[0]^2*
       Co[0]/Ωk[r]^5)^(1/6),
   T[0] == (1/2)*(Ωk[r]*h[0])^2*(μ/
       Rg)*(1/(1 + β[0])),
   Kkr[0] == (k0*(Σ[0]*h[0]))/T[0]^(7/2)};

eq = {Σ'[
     t] == -Σ[
       t] + κ^(-1/5) α^(-4/5) μ^(3/
         5) Ωk[r]^(2/5)*Co[t]^(-1/5),
   h'[t] == -h[
       t] + (κ α Σ[t]^2 Ωk[
          r]^(-5) Co[t])^(1/6),
   T'[t] == -T[t] + 
     1/2 μ/Rg (Ωk[r]^2 h[t]^2)/(1 + β[t]),
   Kkr'[t] == -Kkr[t] + k0 Σ[t]/h[t]*T[t]^(-7/2),
   β'[
     t] == -β[t] + μ/
       Rg (4 σ)/(3 c) T[t]^3/Σ[t] h[t],
   Co'[t] == -Co[t] + (1 + β[t])^4*(1 + Kkr[t]/ke)};



t0 = 20; sol = 
 ParametricNDSolveValue[{eq, ic}, {Σ[t], h[t], 
   T[t]}, {t, 0, t0}, {r}]

Table[Plot3D[sol[r][[i]], {t, 0, t0}, {r, 10^4, 10^5}, 
  PlotRange -> All], {i, 3}]

Table[Plot[sol[r][[i]] /. t -> t0, {r, 10^4, 10^5}, 
  PlotRange -> All], {i, 3}]

fig1

$\endgroup$
  • $\begingroup$ What does t do'? If I change it, the shape of the functions change. $\endgroup$ – mattiav27 Aug 21 at 7:04
  • $\begingroup$ @mattiav27 t is an iterative parameter. Here we use a continuous parameter in the model, but when numerically integrating a system of equations, a discrete parameter is automatically used. $\endgroup$ – Alex Trounev Aug 21 at 9:51
0
$\begingroup$

Not a full answer, just showing how to cast it as a fixed point problem. You have a circular dependency on Co and β, so write this function...

funk[{β_, Co_, r_}] := Module[{Σ, h, T, Kkr, βout, Coout},

  Σ = κ^(-1/5) α^(-4/5) μ^(3/5) Ωk[r]^(2/5)*Co^(-1/5);
  h = (κ α Σ^2 Ωk[r]^(-5) Co)^(1/6);
  T = 1/2 μ/Rg (Ωk[r]^2 h^2)/(1 + β);
  Kkr = k0 Σ/h*T^(-7/2);
  βout = μ/Rg (4 σ)/(3 c) T^3/Σ  h;
  Coout = (1 + βout)^4*(1 + Kkr/ke);

  {βout, Coout, r}
  ];

Now you can run FixedPoint[] or FixedPointList[] and see if it converges. Let it run for 1000 iterations and plot Co. It ends up oscillating between about $10^{54}$ and $0$. The second plot is with the lines through the points.

res = FixedPointList[funk, {.1, 1, 10^4}, 1000];

ListPlot[res[[All, 2]]]
ListLinePlot[res[[All, 2]]]

enter image description here

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.