Trying to use output of ParametricNDSolve as Dirichlet conditions for another ParametricNDSolve

Question

My apologies if the problem was already solved in a similar thread, but i was not able to find a working answer.

I am trying to numerically solve a quite cumbersome PDE of second order which requires numerical initial condition, that I am trying to pass to the NDSolve as Dirichlet condition.

My final goal is to find a function f[x,y] which depends on 2 parameter, and I need to impose initial condition of the form: f[x0,y] = g[y], f[x,y0] = h[x].

Where bot g[y] and h[x] are numerical expressions output of ParametricNDSolve.

Here is the code I am using:

To write the PDE in a simpler form I defined:

a[t_] := t^(2/3);
H[t_] := a'[t]/a[t];
\[Mu][t_, r_] := M/(2*a[t]*r);
DV[t_, r_] := 2*m2*f[t, r];

Then I construct the initial condition:

\[Phi]BG = 
  ParametricNDSolve[{\[Phi]''[t] + 2/t*\[Phi]'[t] + m2*\[Phi][t] == 
     0, \[Phi][1] == 10^{-5}, \[Phi]'[1] == 10^{-5}}, \[Phi], {t, 1, 
    50}, {M, m2}];

\[Psi]Sch = 
  ParametricNDSolve[{ \[Psi]''[r]/(
      1 - M/r) + \[Psi]'[r] (M/r^2 + 2/r (1 - M/r)) + 
      2 m2*\[Psi][r] == 0, \[Psi][20] == 10^{-5}, \[Psi]'[20] == 
     10^{-5}}, \[Psi], {r, 2, 20}, {M, m2}];

\[Phi]r[r_] := 1;
\[Psi]t[t_] := 1;

\[Phi]tot[t_, r_, M_, 
   m2_] := (Evaluate[\[Phi][M, m2][t]] /. \[Phi]BG)*\[Phi]r[r];
\[Psi]tot[t_, r_, M_, 
   m2_] := (Evaluate[\[Psi][M, m2][r]] /. \[Psi]Sch)*\[Psi]t[t];

The last four line are just to define a final 2-D initial condition starting from the 1-D one.

\[Psi]tot[2, 4, 1, 1]
\[Phi]tot[2, 4, 1, 1]

(the above two lines are just to verify that the numerical expression are evaluated correctly).

Finally I wrote down my PDE and ask to evaluate f[1, 1][2, 8]:

sol = ParametricNDSolve[{-((
        (1 + \[Mu][t, r])^2)/(1 - \[Mu][t, r])^2) (D[
          f[t, r], {t, 2}] + 3 H[t]*D[f[t, r], t])  + 
      2 \[Mu][t, 
        r] (4 - 3 \[Mu][t, r])*(1 + \[Mu][t, 
          r])/(1 - \[Mu][t, r])^3 H[t]*D[f[t, r], t] + 
      1/(a[t]^2 (1 + \[Mu][t, r])^4) (D[f[t, r], {r, 2}] + 
         2/(r (1 - \[Mu][t, r]^2)) D[f[t, r], r]) + DV[t, r] == 0 , 
    DirichletCondition[f[t, r] == \[Phi]tot[t, r, M, m2], r == 50] , 
    DirichletCondition[f[t, r] == \[Psi]tot[t, r, M, m2], t == 1]}, 
   f, {t, 1, 10}, {r, 2, 50}, {M, m2}];
    Evaluate[f[1, 1][2, 8]] /. sol

But I obtain the following output:

During evaluation of In[1]:= InterpolatingFunction::dmval: Input value {20.383} lies outside the range of data in the interpolating function. Extrapolation will be used.

During evaluation of In[1]:= InterpolatingFunction::dmval: Input value {21.4043} lies outside the range of data in the interpolating function. Extrapolation will be used.

During evaluation of In[1]:= InterpolatingFunction::dmval: Input value {22.4255} lies outside the range of data in the interpolating function. Extrapolation will be used.

During evaluation of In[1]:= General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation.

During evaluation of In[1]:= SparseArray::drnk: The requested dimensions, {113,1}, have length inconsistent with the tensor rank (3) of the input.

During evaluation of In[1]:= LinearSolve::exopt1: The option setting Method -> Multifrontal cannot be used with arbitrary-precision or exact arguments.

During evaluation of In[1]:= ParametricNDSolve::fempsf: PDESolve could not find a solution.

Out[14]= ParametricFunction[ <> ][1, 1][2, 8]

Any Idea of why I am not obtaining numerical values for f? My guess is that I am not passing correctly the Dirichletconditions, but I don't know how to correct. I also tried something like:

DirichletCondition[f[t, 50] == \[Phi]tot[t, r, M, m2], True] , 
DirichletCondition[f[1, r] == \[Psi]tot[t, r, M, m2], True]},

But did not work either. Am I missing something trivial? Is there any problem in passing initial conditions as output of ParametricNDsolve withouth specifying the parameters before?

Answer 1

There are no sufficient number of boundary and initial conditions for PDE, and boundary and initial conditions are inconsistent. If you resolve this issue then the working code should be like this:

a[t_] := t^(2/3);
H[t_] := a'[t]/a[t];
\[Mu][t_, r_] := M/(2*a[t]*r);
DV[t_, r_] := 2*m2*f[t, r];
\[Phi]BG = 
  ParametricNDSolveValue[{\[Phi]''[t] + 2/t*\[Phi]'[t] + 
      m2*\[Phi][t] == 0, \[Phi][1] == 10^-5, \[Phi]'[1] == 
     10^-5}, \[Phi], {t, 1, 10}, {M, m2}];

\[Psi]Sch = 
  ParametricNDSolveValue[{\[Psi]''[r]/(1 - M/r) + \[Psi]'[
        r] (M/r^2 + 2/r (1 - M/r)) + 2 m2*\[Psi][r] == 
     0, \[Psi][20] == 10^{-5}, \[Psi]'[20] == 10^{-5}}, \[Psi], {r, 2,
     20}, {M, m2}];

\[Phi]r[r_] := 1;
\[Psi]t[t_] := 1;

\[Phi]tot[t_, r_, M_, m2_] := \[Phi]BG[M, m2][t]*\[Phi]r[r];
\[Psi]tot[t_, r_, M_, m2_] := \[Psi]Sch[M, m2][r]*\[Psi]t[t];

sol = ParametricNDSolveValue[{-((0 (1 + \[Mu][t, r])^2)/(1 - \[Mu][t, 
             r])^2) (D[f[t, r], {t, 2}] + 3 H[t]*D[f[t, r], t]) + 
     2 \[Mu][t, 
       r] (4 - 3 \[Mu][t, r])*(1 + \[Mu][t, r])/(1 - \[Mu][t, r])^3 H[
       t]*D[f[t, r], t] + 
     1/(a[t]^2 (1 + \[Mu][t, r])^4) (D[f[t, r], {r, 2}] + 
        2/(r (1 - \[Mu][t, r]^2)) D[f[t, r], r]) + DV[t, r] == 0, 
   f[t, 50] == First[\[Phi]tot[t, 50, M, m2]], 
   f[1, r] == First[\[Psi]tot[1, r, M, m2]]}, 
  f, {t, 1, 10}, {r, 2, 20}, {M, m2}]

Your test Evaluate[f[1,1][2,8]]/.sol now is

sol[1, 1][2, 8]

(*Out[]= 0.00024458*)

Visualisation

Plot3D[sol[1, 1][t, r], {t, 1, 10}, {r, 2, 20}, 
 ColorFunction -> "Rainbow", Mesh -> None, PlotRange -> All, 
 AxesLabel -> Automatic]