use of parametric ndsolve with a constraint

Question

I want to solve differential equation $$ \frac{y''[x]}{(1+y'[x]y'[x])^{3/2}} = -a -y[x]/ \sqrt{2} + x/ \sqrt{2} $$ subject to boundary condition $y(-1) = y(1) = 0$ for some value of $a$. $a$ is found from the condition $\int^{1}_{-1} y[x]\mathrm{d}x = M $ for a given $M$.

This equation arises when one needs to find the shape of a drop resting on a solid surface under gravity.

To start with numerical solution with NDSolve for a fixed value $a$ works fine.

   asymSolNumGrav = 
     NDSolve[{yG''[x]/((1 + yG'[x] yG'[x]))^(3/2) == -2921/10000 - 
         yG[x]/Sqrt[2] + x /Sqrt[2], yG[-1] == 1/10000, yG[1] == 1/1000}, 
      yG[x], {x, -1, 1}, Method -> {"StiffnessSwitching"}, 
      AccuracyGoal -> 20]
    Plot[yG[x] /. asymSolNumGrav, {x, -1, 1}, AspectRatio -> Automatic]
    NIntegrate[yG[x] /. asymSolNumGrav[[1]], {x, -1, 1}]

this gives

Next, I tried to find $a$ given the integral constraint with parametric NDSolve as follows along the lines given in another thread Constraint of NDSolve with an integral of the solution

paraNDsol = 
     ParametricNDSolveValue[{yP''[x]/((1 + yP'[x] yP'[x]))^(3/2) == -a - 
         yP[x]/Sqrt[2] + x /Sqrt[2], yP[-1] == 1/1000, yP[1] == 1/1000}, 
      yP, {x, -1, 1}, {a}, Method -> {"StiffnessSwitching"}, 
      AccuracyGoal -> 20]
M[a_?NumericQ] := (NIntegrate[paraNDsol[a][x], {x, -1, 1}] - 0.2)^2
sol1 = NMinimize[{M[a], 1/10 < a < 99/100}, a]
{Plot[Evaluate[Table[paraNDsol[a][x], {a, 0.1, 1, .1}]], {x, -1, 1}, 
  AspectRatio -> Automatic], 
 Plot[paraNDsol[a][x] /. Last[sol1], {x, -1, 1}, 
  AspectRatio -> Automatic]}

This gives error

 ParametricNDSolveValue::ndsz
    NIntegrate::inumr

How can I get a smooth numerical solution with ParametricNDSolveValue in this case?

Michael Seifert · Accepted Answer · 2021-09-27T16:13:05.770

It appears that Mathematica cannot integrate the ODE all the way from -1 to 1 for certain values of $a$. The important error messages that are returned are

ParametricNDSolveValue::ndsz: At x$8411 == -0.760075, step size is effectively zero; singularity or stiff system suspected.

NMinimize::nnum: The function value (-0.2+NIntegrate[paraNDsol[0.747922][x],{x,-1,1}])^2 is not a number at {a} = {0.747922}.

What's happening here NIntegrate is trying to sample M at various points, including for $a = 0.747922$. It calls ParametricNDSolve, which tries to solve the ODE with this $a$ value and fails, generating the first error and returning a non-numerical answer. NIntegrate then complains that it didn't get a numerical answer.

Thankfully, one can still look at the behavior of M and figure out where the minimum actually is:

Plot[M[a], {a, 0.1, 0.99}, PlotPoints -> 10, MaxRecursion -> 1, PlotRange -> All]

(The options are there to make this code run in a reasonable amount of time.) It appears the function reaches zero between 0.1 and 0.3, so we can search just in this range and avoid the trouble spots:

sol1 = NMinimize[{M[a], 1/10 < a < 3/10}, a]
(* {1.6938110^-21, {a -> 0.19603}} )
Plot[paraNDsol[a][x] /. sol1[[2]], {x, -1, 1}]

What's a bit troubling about this example is that ParametricNDSolve will happily return solutions for which the boundary conditions are not satisfied, without giving any warnings (MM 12.3.1, Mac OSX ARM):

Plot[paraNDsol[0.9][x], {x, -1, 1}]

I would therefore advise caution in using this "graph" technique without actually plotting out your results and making sure that they're doing what you want them to.

use of parametric ndsolve with a constraint

1 Answers1