# Using NeumannValue in finite element method to replace a derivative boundary condition in NDSolve to overcome a computational overflow

I am trying to implement and evaluate the function GeneratePsiOfX. The definition of GeneratePsiOfX is:

GeneratePsiOfX[psi0_?NumericQ,xmin_?NumericQ,xmax_?NumericQ,maxsteps_?NumericQ,method_]:=Module[{ode,ID,ORD,eq1,ss,as,Eqs,Soln,Soln2,tmp},
ode=2 (\[Psi]')[x]+x (\[Psi]'')[x]+3x E^\[Psi][x] Sqrt[\[Pi]] Erf[Sqrt[\[Psi][x]]]-2 Sqrt[\[Psi][x]] (3+2 \[Psi][x])x;
ID={\[Psi]0->psi0,\[Psi]p0->0};
ORD=10;
eq1=\[Psi]->Function[x,Sum[a[i]x^i,{i,0,ORD}]/.{a[0]->\[Psi]0,a[1]->\[Psi]p0}/.ID];
ss=Series[ode/.eq1,{x,0,ORD}];
(*Print[ss];*)
as=Table[a[i],{i,2,ORD}];
Eqs=Table[SeriesCoefficient[ss,i]==0,{i,1,ORD}];
(*Print[Eqs];*)
seriesSoln=Solve[Eqs,as]//Flatten;
Soln=\[Psi][x]/.eq1/.seriesSoln//N;
ode1=2 (\[Psi]')[x]+x (\[Psi]'')[x]+3x E^\[Psi][x] Sqrt[\[Pi]] Erf[Sqrt[\[Psi][x]]]-2 Sqrt[\[Psi][x]] (3+2 \[Psi][x])x==0;
xstart=xmin;
xstop=xmax;
Soln2=NDSolve[{ode1,\[Psi][xstart]==Soln/.x->xstart,\[Psi]'[xstart]==D[Soln,x]/.x->xstart},\[Psi][x],{x,xstart,xstop},MaxSteps->maxsteps,Method->method][[1,1,2]];
Piecewise[{{Soln,x<=xstart},{Soln2,x>xstart}}]
]


Upon evaluating GeneratePsiOfX as follows:

GeneratePsiOfX[15,0.001,4,10^7,{"StiffnessSwitching"}]


I get the following error:

General::ovfl: Overflow occurred in computation.


So, I decided to switch to Finite Element Method instead of Stiffness Switching as below:

GeneratePsiOfX[15,0.001,4,10^7,{"FiniteElement"}]


I now get the following error:

NDSolve::fembderiv: The expression (\[Psi]')[0.001]==-29568.3 given as a spatial boundary condition for the possibly automatically chosen finite element method should not have explicit derivatives of the dependent variables. NeumannValue should be used to specify spatial derivatives at the boundary.


I read the official documentation of NeumannValue but it is a little difficult to understand. I tried reading a few posts on this website as well. Here is what I think is required. I modified the NDSolve function within the GeneratePsiOfX function to include NeumannValue as follows:

Soln2 = NDSolve[{ode1, \[Psi][xstart] == Soln /. x -> xstart, NeumannValue[D[Soln,x]/.x->xstart,x==xstart]}, \[Psi][x], {x, xstart, xstop}, MaxSteps -> maxsteps,
Method -> method][[1, 1, 2]];


However, when I implement the following:

GeneratePsiOfX[15, 0.001, 4, 10^7, {"FiniteElement"}]


I get the error:

NDSolve::deqn: Equation or list of equations expected instead of NeumannValue[-29568.3,x==0.001] in the first argument {3 E^\[Psi][x] Sqrt[\[Pi]] x Erf[Sqrt[\[Psi][x]]]-2 x Sqrt[\[Psi][x]] (3+2 \[Psi][x])+2 (\[Psi]')[x]+x (\[Psi]'')[x]==0,\[Psi][0.001]==10.8068,NeumannValue[-29568.3,x==0.001]}.


So, Mathematica first tells me to use NeumannValue but then when I do so, it tells me that's not what it expects. I am a little confused as what to do. Also, I am using Mathematica v.13.0.1.0.

• As the error message suggests, NeumannValue should be used as part of an equation, and not on its own as you have it. Jul 3, 2022 at 23:18
• In addition to @MarcoB's comment: 1. Most of all, I don't think Overflow is something can be overcomed by method changing, better to doubt-check if the equation itself is correct or not. 2. FiniteElement method isn't designed for IVP, if you insist on using it, check this post: mathematica.stackexchange.com/q/172972/1871 Jul 4, 2022 at 2:11
• @MarcoB Would you be able to suggest how to use NeumannValue as part of the ODE equation. I tried ic2 = NeumannValue[D[Soln, x] /. x -> xstart, x == xstart]; and ode1=2 (\[Psi]^\[Prime])[x]+x (\[Psi]^\[Prime]\[Prime])[x]+3x E^\[Psi][x] Sqrt[\[Pi]] Erf[Sqrt[\[Psi][x]]]-2 Sqrt[\[Psi][x]] (3+2 \[Psi][x])x == ic2; and Soln2 = NDSolve[{ode1,\[Psi][xstart]==Soln/.x->xstart},\[Psi][x],{x,xstart,xstop},MaxSteps->maxsteps,Method->method][[1,1,2]];. But that didn't work. Jul 4, 2022 at 15:14
• @xzczd Overflow doesn't occur for values smaller than 15 for the first argument (psi0) of the function GeneratePsiOfX. So, I think that the ODE equation may be fine. Also, I can' locate it at the moment but I saw a post in this site that suggested using FiniteElement to overcome the Overflow problem. Jul 4, 2022 at 15:56
• @xzczd Here is the post where user21 uses FiniteElement method. The question was regarding an Overflow issue. I hope I interpreted this correctly. Jul 4, 2022 at 16:46