2
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ why are you multiplying the first part by zero there? 0 (1 + \[Mu][t, r])^2)/(1 - \[Mu][t, r])^2) as anything times zero is zero. $\endgroup$
    – Nasser
    Apr 18, 2020 at 20:48
  • $\begingroup$ Thanks for the ready anwer, that actually is actually a typo from copying, was not present in the original code. I will edit it, thanks againt $\endgroup$
    – Leonardo
    Apr 18, 2020 at 21:45

1 Answer 1

4
$\begingroup$

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]

Figure 1

$\endgroup$
4
  • $\begingroup$ I like your answer more than mine. So deleting mine. But note that still the M and m2 parameters range is needed to be known. For example Evaluate[sol[4, 2][1, 1]] gives 1/0 error. But I think your solution looks better than what I show. $\endgroup$
    – Nasser
    Apr 18, 2020 at 21:29
  • $\begingroup$ @Nasser Sorry, I didn't read your post attentively, only check your notes about 10^{-5} and made some corrections in my code. Thank you. Also it is better to use ParametricNDSolveValue for this kind of problem. $\endgroup$ Apr 18, 2020 at 21:40
  • $\begingroup$ Thanks very much for the answer, it was very helpful. If I understand properly there was no need of use DirichletCondition, but for sake of clarity, it was the use of it the problem or the specific choice initial condition? $\endgroup$
    – Leonardo
    Apr 18, 2020 at 21:54
  • 1
    $\begingroup$ @Leonardo If we use DirichletCondition then it means that we force NDSolve to use FEM as a method of solution. But FEM is a very sensitive and not applicable to some problems including this one. $\endgroup$ Apr 18, 2020 at 22:03

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.