3
$\begingroup$

I am trying to solve the system of equations:

eq1 = (Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 - 
   2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[λ[t, r], r])/
    R[t, r] + 
   2 (Exp[-2 λ[t, r]] D[R[t, r], {r, 2}])/R[t, r] - (
   Exp[-2 ν[t, r]] D[R[t, r], t]^2)/R[t, r] - 
   2 (Exp[-2 ν[t, r]] D[R[t, r], t] D[ν[t, r], t])/R[t, r] + (
    8 Pi G)/3 ϵ[t, r];

eq3 = (Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 + 
   2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[ν[t, r], r])/
    R[t, r] - (Exp[-2 ν[t, r]] D[R[t, r], t]^2)/R[t, r]^2 + 
   2 (Exp[-2 ν[t, r]] D[R[t, r], t] D[ν[t, r], t])/R[t, r] - 
   2 (Exp[-2 ν[t, r]] D[R[t, r], {t, 2}])/R[t, r] - (8 Pi G)/
    3 ϵ[t, r]/3;

where the functions $\epsilon$ and $\nu$ are given by:

ϵ[t_, r_]=(E^( - 4/3 λ[t, r]) (1/100)^(8/3))/R[t, r]^(8/3)
ν[t_, r_]=1/4 (Log[1/(1 + t)^(8/3)] - Log[E^(-(4/3) λ[t, r])/R[t, r]^(8/3)])

I use the method pdetoode by @xzczd as follows:

lb = 10^-1;
rb = 1;
xdifforder1 = 2;
xdifforder2 = 1;
ic1 = {R[0, r] == 10^-2, Derivative[1, 0][R][0, r] == 1};
ic2 = {λ[0, r] == 0};
bc1 = {R[t, 0] == t + 10^-2, R[t, 10^-1] == 10^-2, 
   Derivative[0, 1][R][t, 0] == 0};
bc2 = {λ[t, 0] == 0};
points = 25;
grid = Array[# &, points, {lb, rb}];
removeredundant1 = #[[2 ;; -2]] &;
removeredundant2 = #[[1 ;; -1]] &;
ptoofunc1 = pdetoode[R[t, r], t, grid, xdifforder1];
ptoofunc2 = pdetoode[λ[t, r], t, grid, xdifforder2];
odeqn1 = eq1 // ptoofunc1 // removeredundant1;
odeqn2 = eq3 // ptoofunc2 // removeredundant2;
odeic1 = removeredundant1 /@ ptoofunc1@ic1;
odeic2 = removeredundant2 /@ ptoofunc2@ic2;
odebc1 = bc1 // ptoofunc1;
odebc2 = bc2 // ptoofunc2;
sollst = NDSolveValue[{odebc1, odebc2, odeic1, odeic2, odeqn1, odeqn2}, {R /@ grid, [Lambda] /@ grid}, {t, 0, 1},  MaxSteps -> Infinity]
sol = rebuild[sollst, grid]

but I get the error:

enter image description here

So it seems that while it works well for $R$, the grid points are not mapped onto $\lambda$ and it becomes a function of a list. I do not understand why it happens, since it seems to me that I am applying the same method to both $R$ and $\lambda$. Can anyone help?

EDIT

After the comment by @xzczd, my code now looks like this:

lb = 10^-1;
rb = 1;
xdifforder = 4;
ic = {R[0, r] == 10^-2, 
   Derivative[1, 0][R][0, r] == 1, λ[0, r] == 0};
bc = {R[t, 0] == t + 10^-2, R[t, 10^-1] == 10^-2, 
   Derivative[0, 1][R][t, 0] == 0, λ[t, 0] == 0};
points = 25;
grid = Array[# &, points, {lb, rb}];
removeredundant = #[[3 ;; -3]] &;
ptoofunc = pdetoode[{R[t, r], λ[t, r]}, t, grid, xdifforder];
odeqn1 = eq1 // ptoofunc // removeredundant;
odeqn2 = eq3 // ptoofunc // removeredundant;
odeic = removeredundant /@ ptoofunc@ic;
odebc = bc // ptoofunc;

but now Mathematica complaints that there are fewer variables than equations

$\endgroup$
16
  • $\begingroup$ Why do you use pdetoode twice? $\endgroup$
    – xzczd
    Commented May 20, 2021 at 7:02
  • $\begingroup$ @xzczd because I have two functions $R$ and $\lambda$, I have found some posts on this site for system of equations and, as far as I understand, this seems to be the method. Is it wrong? $\endgroup$
    – mattiav27
    Commented May 20, 2021 at 7:05
  • $\begingroup$ Yes, it's wrong. If the domain of the unknown functions are the same, then you only need to use pdetoode once, i.e. pdetoode[{R[t, r], \[Lambda][t, r]}, … or simply pdetoode[{R, \[Lambda]}[t, r], … $\endgroup$
    – xzczd
    Commented May 20, 2021 at 7:10
  • $\begingroup$ @xzczd but the functions have different differentiation order: $R$ it is 2, while $\lambda$ it is 1 $\endgroup$
    – mattiav27
    Commented May 20, 2021 at 7:13
  • $\begingroup$ diff is short for difference. You may want to read the comments here: mathematica.stackexchange.com/a/174775/1871 $\endgroup$
    – xzczd
    Commented May 20, 2021 at 7:21

1 Answer 1

3
$\begingroup$

This system can be solved with DAE solver on the short time interval {t,0,.08} as follows

ϵ = (Exp[(-4/3 λ[t, r])] (1/100)^(8/3))/R[t, r]^(8/3);
ν = 1/
    4 (Log[1/(1 + t)^(8/3)] - 
     Log[E^(-(4/3) λ[t, r])/R[t, r]^(8/3)]);

eq = {(Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 - 
    2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[λ[t, r], r])/
      R[t, r] + 
    2 (Exp[-2 λ[t, r]] D[R[t, r], {r, 2}])/
      R[t, r] - (Exp[-2 ν] D[R[t, r], t]^2)/R[t, r] - 
    2 (Exp[-2 ν] D[R[t, r], t] D[ν, t])/R[t, r] + (8 Pi G)/
      3 ϵ,
   (Exp[-2 λ[t, r]] D[R[t, r], r]^2)/R[t, r]^2 + 
    2 (Exp[-2 λ[t, r]] D[R[t, r], r] D[ν, r])/
      R[t, r] - (Exp[-2 ν] D[R[t, r], t]^2)/R[t, r]^2 + 
    2 (Exp[-2 ν] D[R[t, r], t] D[ν, t])/R[t, r] - 
    2 (Exp[-2 ν] D[R[t, r], {t, 2}])/R[t, r] - (8 Pi G)/
      3 ϵ/3} // Simplify

ic = {R[t, r] == 10^-2, D[R[t, r], t] == 1, λ[t, r] == 0} /. 
   t -> 0;
bc = {R[t, r] - t - 10^-2 == 0, 
    D[R[t, r], r] == 0, λ[t, r] == 0} /. r -> lb;


G = 1; lb = 10^-1;
rb = 1;
sol = NDSolve[{eq[[1]]+eq[[2]] == 0, eq[[2]] == 0, ic, 
   bc}, {R, λ}, {t, 0, .08}, {r, lb, rb}, 
  Method -> {"IndexReduction" -> Automatic, 
    "EquationSimplification" -> "Residual", 
    "PDEDiscretization" -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MinPoints" -> 80, "MaxPoints" -> 80, 
        "DifferenceOrder" -> 8}}}]

{Plot3D[Evaluate[R[t, r] /. sol[[1]]], {t, 0, 0.08}, {r, lb, rb}, 
  ColorFunction -> Hue, Mesh -> None, AxesLabel -> Automatic, 
  PlotRange -> All, Boxed -> False, PlotPoints -> 50], 
 Plot3D[Evaluate[λ[t, r] /. sol[[1]]], {t, 0, 0.08}, {r, lb, 
   rb}, ColorFunction -> Hue, Mesh -> None, AxesLabel -> Automatic, 
  Boxed -> False, PlotPoints -> 50]}

It looks very regular, but it is not clear why system blow up after t=0.089 Figure 2

There is a stable solution for $\lambda >0$ with the set of boundary conditions

ic = {R[t, r] == 10^-2, 
    D[R[t, r], t] == 1/100, λ[t, r] == 0} /. t -> 0;
bc = {R[t, r] - t/100 - 10^-2 == 0, 
    D[R[t, r], r] == 0, λ[t, r] == 0} /. r -> lb;

This solution on {t,0,1} looks like

Figure 2

Note, we are not been able to reproduce any of these results with using pdetoode or pdetoae.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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