# Problem with a solution of PDE with initial and boundary conditions

I would like to solve relativistic hydrodynamic equations (nonlinear PDEs) introduced here:

I use eqs (33 - 35), (38 - 41), where (40) P(rho)=k*rho^g0 (all with one spatial coordinate "mu" and one temporal "t"). CODE EDITED: 16.04.2020 - this I use.

 (*Initial functions-stationary,homogeneous perfect fluid sphere \
structure*)
(****************************************************************)

ClearAll["Global*"]
Needs["NDSolveFEM"]
c = 2.99792*10^10;(*m/s*)
gr = 6.674*10^-8;(*grav. const. in cm^3*g^-1*s^-2*)
gcc = gr/c^2;
m0 = 1.672621*10^-24*gr/c^2;(*proton mass in g trnasformed to cm*)
Ms0 = 1.98855*10^33;
Ms = Ms0*gr/c^2;(*mass of central object in g trnasfomred to cm*)
dr = 10^-5;(*small step and initial m is only e*)
(*initital data*)
g0 = 5/3; rho0 = 10^11; ep0 = 3.64*10^18; e0 =
rho0 (1 + ep0/c^2); pc = (g0 - 1)*rho0*ep0;
dmu = 4*\[Pi]*rho0*dr^2; mumax = 21 Ms0; \[Gamma] = g0; k = pc/rho0^g0;
{pc // N, rho0 // N, e0, ep0 // N, ep0/c^2}

(*Solution TOV and mass equation*)
s = NDSolve[{r'[mu] == Sqrt[1 - 2 m[mu]*gr/(r[mu]*c^2)]/(
4 \[Pi]*rho0*r[mu]^2),
m'[mu] == e0/rho0 Sqrt[1 - (2 m[mu] gcc)/r[mu]], r[dmu] == dr,
m[dmu] == dmu}, {r, m}, {mu, dmu, mumax}];
(*Initial functions to hydrodynamical calculations*)
r0 = r /. s[[1, 1]]; fm0 = m /. s[[1, 2]];
{r0[mumax], fm0[mumax]/Ms0, dmu // N, mumax // N}
f3 = Plot[{fm0[mu]}/Ms0, {mu, dmu, mumax}, Frame -> True,
FrameLabel -> {"\[Mu] [g]", "M/Ms []"}, PlotRange -> All]
f4 = Show[
Plot[{r0[mu]}, {mu, dmu, mumax}, Frame -> True,
FrameLabel -> {"\[Mu] [g]", "r [cm]"}]]
frho0[x_] = If[x < mumax, rho0, 1];

(*Relativistic hydrodynamical equations-collapse of star*)
(**************************************************)

(*introducing of equation*)
G[mu_, t_] = 4 \[Pi]*rho[mu, t]*r[mu, t]^2*D[r[mu, t], mu];(*MW39*)
w[mu_, t_] = 1 + ep[mu, t]/c^2 + p[mu, t]/(rho[mu, t]*c^2);(*MW41*)
a[mu_, t_] = 1/w[mu, t];
ep[mu_, t_] = k*rho[mu, t]^(\[Gamma] - 1)/(\[Gamma] - 1);
p[mu_, t_] = (\[Gamma] - 1) ep[mu, t]*rho[mu, t];(*MW40*)
equt[mu_,
t_] = -a[mu,
t] (4 \[Pi]*r[mu, t]^2*G[mu, t]/w[mu, t]*D[p[mu, t], mu] + (
m[mu, t]*gr)/
r[mu, t]^2 + (4 \[Pi]*gr)/c^2 p[mu, t]*r[mu, t]);(*MW33*)
eqrt[mu_, t_] = a[mu, t]*u[mu, t];(*MW34*)
eqmm[mu_, t_] =
4 \[Pi]*rho[mu, t]*(1 + ep[mu, t]/c^2)*
r[mu, t]^2 D[r[mu, t], mu];(*MW38*)
eqrhort[mu_, t_] = -a[mu, t]*rho[mu, t]*r[mu, t]^2 D[u[mu, t], mu]/
D[r[mu, t], mu];(*MW35*)

(*preparation for solution*)
(*equations*)
eqs = {D[u[mu, t], t] == equt[mu, t], D[r[mu, t], t] == eqrt[mu, t],
D[m[mu, t], mu] == eqmm[mu, t],
D[rho[mu, t]*r[mu, t]^2, t] == eqrhort[mu, t]};
(*boundary conditions*)
bcon = {DirichletCondition[u[mu, t] == 0., mu == dmu],
DirichletCondition[r[mu, t] == r0[dmu], mu == dmu],
DirichletCondition[m[mu, t] == fm0[dmu], mu == dmu],
DirichletCondition[rho[mu, t] == frho0[mumax], mu == mumax]};
(*initial conditions*)
incon = {u[mu, 0] == 0., r[mu, 0] == r0[mu], m[mu, 0] == fm0[mu],
rho[mu, 0] == frho0[mu]};

(*PDEs solution*)
Clear[fu, fr, fm, fro]
{fu, fr, fm, fro} =
NDSolveValue[{eqs, incon, bcon}, {u, r, m, rho}, {mu, dmu,
mumax}, {t, 0, 0.1}]


Initial functions r0[mu], fm0[mu] and frho0[mu] are interpolated functions coming from numerical solution of stationary problem. Result of this solution are error messages:

NDSolveValue::femcnsd: The PDE coefficient -((6.674*10^-8 m[mu])/r[mu]^2)-1.15712*10^-17 r[mu] rho[mu]^(5/3)-3.26355*10^23 r[mu]^4 rho[mu]^(2/3) (r^\[Prime])[mu] (rho^\[Prime])[mu] does not evaluate to a numeric scalar at the coordinate {2.08798*10^34}; it evaluated to Indeterminate instead.
NDSolveValue::femcnsd: The PDE coefficient -((6.674*10^-8 m[mu])/r[mu]^2)-1.15712*10^-17 r[mu] rho[mu]^(5/3)-3.26355*10^23 r[mu]^4 rho[mu]^(2/3) (r^\[Prime])[mu] (rho^\[Prime])[mu] does not evaluate to a numeric scalar at the coordinate {2.08798*10^34}; it evaluated to Indeterminate instead.


Unfortunately, I don't know where is the problem (whole concept, method or...). The problem appears ever in half value of endpoint of integration (mumax/2), doesn't matter what "mumax" is. I'm able to draw (and evaluate in all point of the range) all defined functions in the initial time without problems.

Thank you for help or suggestions.

PS: I'm new here if something is misspelt, marked or unlisted. Please notify me. Thank you.

• I think you will need to provide the missing parameters ($k$, $gr$, $c$, $dmu$, $mumax$, etc) to receive help. Apr 15, 2020 at 13:34
• @Vrbic It looks like gas dynamic equations. It could be better to explain the problem you want to solve. Apr 15, 2020 at 14:39
• @Vrbic it looks like $dr$ is not assigned in dmu = 4*\[Pi]*rho0*dr^2. Apr 15, 2020 at 16:03
• @Tim Laska my mistake when I copied eqs, dr = 10^-5; Apr 15, 2020 at 16:10
• @Vrbic Unfortunately, I think you will need to provide a Minimal Working Example to get more help. One way to check is to start a new Mathematica session and copy the code provided here and see if you obtain the same claimed output. If not, there is more work to do. We can't get very far without your previously calculated interpolation functions. Also, I believe you will need to consider scaling your variables since your $mu$ independent variable ranges over 30 orders of magnitude. Apr 15, 2020 at 16:56

## 1 Answer

First part of code can be used as it is with small modification only. But the last part we should rebuild from the ground. Thanks to paper May & White I found some combination of equations to solve it with NDSolve. All variables in this code should be normalised including t and mu as c*t and mu/mumax. This code allow us to solve up to tm=2.9*10^4 (at this moment initial density increased of 120 times).

c = 2.99792*10^10;(*m/s*)gr =
6.674*10^-8;(*grav.const.in cm^3*g^-1*s^-2*)gcc = gr/c^2;
m0 = 1.672621*10^-24*
gr/c^2;(*proton mass in g trnasformed to cm*)Ms0 = 1.98855*10^33;
Ms = Ms0*gr/c^2;(*mass of central object in g trnasfomred to cm*)dr =
10^-5;(*small step and initial m is only e*)(*initital data*)g0 =
5/3; rho0 = 10^11; ep0 = 3.64*10^18; e0 =
rho0 (1 + ep0/c^2); pc = (g0 - 1)*rho0*ep0;
dmu = 4*\[Pi]*rho0*dr^2; mumax = 21 Ms0; \[Gamma] = g0; k = pc/rho0^g0;
{pc // N, rho0 // N, e0, ep0 // N, ep0/c^2}

(*Solution TOV and mass equation*)
{r0, fm0} =
NDSolveValue[{r'[mu] ==
Sqrt[1 - 2 m[mu]*gr/(r[mu]*c^2)]/(4 \[Pi]*rho0*r[mu]^2),
m'[mu] == e0/rho0 Sqrt[1 - (2 m[mu] gcc)/r[mu]], r[dmu] == dr,
m[dmu] == dmu}, {r, m}, {mu, dmu, mumax}];
(*Initial functions to hydrodynamical calculations*)
frho0[x_] = 1 + rho0 (1 - Tanh[10 (x - .9)])/2;
{r0[mumax], fm0[mumax]/Ms0, dmu // N, mumax // N}
{Plot[fm0[mu], {mu, dmu, mumax}, Frame -> True,
FrameLabel -> {"\[Mu] [g]", "M"}, PlotRange -> All],
Plot[r0[mu], {mu, dmu, mumax}, Frame -> True,
FrameLabel -> {"\[Mu] [g]", "r [cm]"}],
Plot[frho0[mu], {mu, 0, 1}, Frame -> True,
FrameLabel -> {"\[Mu] [g]", "rho"}, PlotRange -> All]}


Parameters scales to normalise parameters

{rhoN, rN, mN, eN,uN} = {rho0 // N, r0[mumax], fm0[mumax],
10^-4 c^2,c};


Relativistic hydrodynamical equations - collapse of star

G[mu_, t_] :=
4 \[Pi]*(rhoN rN^3)*rho[mu, t]*r[mu, t]^2*
D[r[mu, t], mu]/mumax(*MW39*);
p[mu_, t_] := (\[Gamma] - 1) (eN rhoN) ep[mu, t]*rho[mu, t](*MW40*);
w[mu_, t_] :=
1 + eN ep[mu, t]/c^2 + p[mu, t]/(rho[mu, t]*rhoN c^2)(*MW41*);
(*introducing of equation*)
eq = {D[u[mu, t],
t] == (-a[mu,
t] (4 \[Pi] rN^2*r[mu, t]^2*G[mu, t]/w[mu, t]*
D[p[mu, t], mu]/mumax + (m[mu, t]*gr mN/rN^2)/
r[mu, t]^2 + (4 \[Pi]*gr rN)/c^2 p[mu, t]*r[mu, t]))/
c^2(*MW33*), D[r[mu, t], t] == a[mu, t]*u[mu, t](*MW34*),
D[rho[mu, t] r[mu, t]^2, t] == -a[mu, t]*rho[mu, t]*
r[mu, t]^2 D[u[mu, t], mu]/D[r[mu, t], mu]/rN(*MW35*),
D[ep[mu, t], t] == -p[mu, t]/(eN rhoN) D[1/rho[mu, t], t](*36*),
D[a[mu, t] w[mu, t], t] ==
a[mu, t] (D[ep[mu, t], t] eN + p[mu, t] D[1/rho[mu, t], t]/rhoN)/
c^2(*MW37t*),
D[m[mu, t], t] == -4 \[Pi]* rN^3 /mN *p[mu, t]*
r[mu, t]^2 D[r[mu, t], t]/c^2(*MW12*)};


Variables, initial and boundary conditions

var = {rho, r, ep, u, a, m};

{dmu1, mumax1} = {dmu, mumax}/mumax;

ic = {u[mu, 0] == 0., r[mu, 0] == r0[mu mumax]/rN,
m[mu, 0] == fm0[mu mumax]/mN, rho[mu, 0] == frho0[mu ]/rhoN,
a[mu, 0] == 1, ep[mu, 0] == 1};
bc = {u[dmu1, t] == 0.0, r[dmu1, t] == r0[dmu]/rN,
m[mumax1, t] == fm0[mumax]/mN,
rho[mumax1, t] == frho0[mumax1]/rhoN, a[mumax1, t] == 1,
ep[mumax1, t] == 1};


Equations solving and visualisation

tm = 2.5 10^4; Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{frho, fr, fep, fu, fa, fm} =
NDSolveValue[{eq, ic, bc}, var, {mu, dmu1, mumax1}, {t, 0., tm},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 101, "MaxPoints" -> 101,
"DifferenceOrder" -> 2}},
EvaluationMonitor :> (currentTime = t;)];]

{DensityPlot[rho0 frho[mu, t], {mu, dmu1, mumax1}, {t, 0., tm},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "rho", AxesLabel -> Automatic, PlotRange -> All],
DensityPlot[rN fr[mu, t], {mu, dmu1, mumax1}, {t, 0., tm},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "r", AxesLabel -> Automatic, PlotRange -> All],
DensityPlot[c fu[mu, t], {mu, dmu1, mumax1}, {t, 0., tm},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "u", AxesLabel -> Automatic, PlotRange -> All],
DensityPlot[ fa[mu, t], {mu, dmu1, mumax1}, {t, 0., tm},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "a", AxesLabel -> Automatic, PlotRange -> All],
DensityPlot[mN fm[mu, t], {mu, dmu1, mumax1}, {t, 0., tm},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "m", AxesLabel -> Automatic, PlotRange -> All]}


• Big thanks, it works! I managed to reformulate the initial state (so that rho was not a step function). I have to this a few questions: 1) Could you briefly comment on why you chose the MethodOfLines method and subsequent options? 2) Scaling of values is necessary? If so then why? 3) What units does time have? Is it that tm = 1 means end time is 1s? 4) Is it possible to control the time step? Or how is the time step chosen during integration? (I am pointing to the Courant conditions - viz paper) 5) Is it possible to set to preserve an invariant for this solution procedure? Apr 22, 2020 at 8:53
• @Vrbic 1. There is practically one method to solve this problem, therefore we use option Method -> {"MethodOfLines", ...} . 2. Scaling is necessarily step, so never use values like 10^34. 3. Scaled time is c t, and unscaled maximal time is tmax=tm/c. For normal star density it could be better to rescale time as c t/mumax. 4. Time step is controlled automatically, but we can control "MaxPoints"`. 5. We could include invariant in the system and check how it save. Apr 22, 2020 at 10:40
• @ Alex Trounev Thank you again! It helped me a lot. I'm going to play with it :-) Apr 22, 2020 at 11:57
• @Vrbic You are welcome! By the way, how quick is a test code on your computer? Apr 22, 2020 at 12:33
• @ Alex Trounev, the first compilation about 0.5s, the second compilaton about 0.33s - i7,16GB, Win10. Apr 22, 2020 at 16:10