Mathematica 12, Windows 10.

I am trying to solve a PDE in one spatial dimension $R$ and time $t$. I need a solution for displacement $r(R,t)$, radial Cauchy stress $T_{RR}(R,t)$, and radial growth $\gamma_{RR}(R,t)$.

To make things as simple as possible, I first write the static version of the system (two coupled ODEs), making sure we get a solution. Then I add a trivial time evolution turning it into a PDE, and hopefully we can figure out the problem together.

Static problem (coupled ODEs)

 γRR[R_, t_] := 3 
γθθ[R_, t_] := 4

eqns = {
  r[R, t] D[r[R, t], R] == γRR[R, t] γθθ[
     R, t] R,
  r[R, t]^4  R γθθ[R, t] D[TRR[R, t], R] == 
   2 γRR[R, 
     t] (r[R, t]^4 - R^4 γθθ[R, t]^4) 

ic = {r[0.1, t] == 0.1, TRR[1, t] == 0}

sol = First@NDSolve[{eqns, ic}, {r[R, t], TRR[R, t]}, {R, 0.1, 1}]

PlotSol[sol_] := 
 Plot[{r[R, t], TRR[R, t]} /. sol /. t -> 1 // Evaluate, {R, 0.1, 1}, 
  AxesLabel -> {"R", None}, 
  PlotLegends -> Placed[{"r(R)", "TRR(R)"}, {Center, Top}]]

The solution is as expected:

enter image description here

Dynamic problem (PDE)

Now I simply add a trivial equation to the system of differential equations: The time derivative of $\gamma_{RR}$ is zero, that is $\dot{\gamma}_{RR}=0$. I amend the initial conditions with $\gamma_{RR}(R, t=0)=3$, which of course completely decouples and should give the same result as in the ODE example. I now NDSolve now over a domain in both $R$ and $t$:


AppendTo[eqns, D[γRR[R, t], t] == 0]

AppendTo[ic, γRR[R, 0] == 3 ]

sol = First@
  NDSolve[{eqns, ic}, {r[R, t], TRR[R, t]}, {R, 0.1, 1}, {t, 0, 1}]

However, unexpectedly, I get some error messages from the NDSolve and the correct solution is not returned:

enter image description here

What is going on here? I added a trivial uncoupled set of equations which should not have changed the results.


1 Answer 1


It's because NDSolve is using different method for these 2 problems. For the first problem, Shooting method, which is for nonlinear boundary value problem of ordinaray differential equation, is used; while for the second problem, FiniteElement method, which supports nonlinear problem as of v12, is used. Different solvers of NDSolve are not equally powerful. For the second problem, we can obtain a result by adjusting the option a little. (The result is not perfect though):

sol = First@
  NDSolve[{eqns, ic}, {r[R, t], TRR[R, t], γRR}, {R, 0.1, 1}, {t, 0, 1},
          InitialSeeding -> {r[R,t] == 3, TRR[R,t] == 3, γRR[R,t] == 2},
          Method -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.05}]

enter image description here

"Then, is it possible to avoid using pure FiniteElementin NDSolve?" Well, for this specific toy problem, if we differentiate the first 2 equations in t direction and add proper initial condition, it's possible:

γθθ[R_, t_] := 4
tend = 1;
bleft = 1/10; bright = 1;
With[{r = r[R, t], TRR = TRR[R, t], γRR = γRR[R, t], γθθ = γθθ[R, t]},

  eqns = {r D[r, R] == γRR γθθ R,
          r^4 R γθθ D[TRR, R] == 2 γRR (r^4 - R^4 γθθ^4),
          D[γRR, t] == 0};

  bc = {r == 0.1 /. R -> bleft, TRR == 0 /. R -> bright};

  ic = γRR == 3 /. t -> 0;];

Block[{γRR}, γRR[R_, t_] := 3; 
  solinitlst = 
   NDSolveValue[{eqns // Most, bc}, {r[R, t], TRR[R, t]}, {R, bleft, bright}]];

icadd = {r[R, 0], TRR[R, 0]} == solinitlst;

sollst = NDSolveValue[{D[eqns // Most, t], eqns // Last, bc, ic, icadd}, {r, 
    TRR, γRR}, {t, 0, tend}, {R, bleft, bright}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 50, 
       "MinPoints" -> 50, "DifferenceOrder" -> 2}}];

Plot[sollst[R, tend] // Through // Evaluate, {R, bleft, bright}]

enter image description here

Notice this approach will probably fail in more complicated cases, because the DAE solver NDSolve is not that strong. In those cases, advanced technique is needed, here is an example.

  • 2
    $\begingroup$ (+1), it's usually better to solve time dependent problems with the method of lines and FEM/TGP as a spatial discretization then as a pure spatial problem like done here. Usually specifying Method->"MethofOfLInes" then indicates some setup problem for the time dependent problem and is ultimately the reason a pure spatial FEM as chosen. $\endgroup$
    – user21
    Dec 19, 2019 at 9:50
  • $\begingroup$ So, for a problem like this, would it be best to use NDSolve only for the spatial problem (say, with default properties, like in my ODE example) and then write a loop (say simple Euler scheme) iterating over time? Or is there a magic NDSolve trick that could specify different methods for spatial and temporal integration as options? $\endgroup$ Dec 26, 2019 at 14:59
  • $\begingroup$ @AlexanderErlich "……Euler scheme……" I doubt, because we need to transform the system to the standard form required by Euler scheme, which is not easy for the specific problem. (Related: mathematica.stackexchange.com/a/158519/1871) “…magic NDSolve trick…" There doesn't seem to be, because NDSolve is having difficulty in noticing MethodOfLines can be used in t direction. (Adding Method -> {MethodOfLines, TemporalVariable -> t} causes tvic warning. ) $\endgroup$
    – xzczd
    Dec 27, 2019 at 4:15
  • 1
    $\begingroup$ @AlexanderErlich Anyway, I've added a solution for the toy example forcing NDSolve to use DAE solver in t direction. It'll probably fail for more complicated problem, though. $\endgroup$
    – xzczd
    Dec 27, 2019 at 5:17

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.