1
$\begingroup$

I am trying to solve a 1D non-linear ODE with FEM. The ODE corresponds to a Model with 3 parameters ([Beta], [lambda], A_2). The ODE reads

mpl/(m^2)Cos(z)^4 f''(z)-2mpl/(m^2)Cos(z)^4Tan(z)f'(z)-A2\rho(z)/mpl *f(z) + Exp(-[lambda]f(z)/mpl+Log[[lambda]+Log[10][Beta]),

where mpl, m are some fixed constants, \rho(z) is a density. I have implemented the above equation and Mathematica does return a solution for the parameters that are of interest to me. I give an initial seed to NDSolve, and create my own mesh.

For many parameters it does return the correct solution (I have a very good analytical approximation to check that) but for some it hardly improves on the not very accurate seed and still returns a solution.

My question is: For the parameters [beta]=1, A2=10^(40), \lambda=10^27; it does return the correct solution, for [beta]=1, A2=10^(40), \lambda=10^(31) it does not and seems to get stuck on the seed. What is the reason for this? I provide the full code necessary to just run the final differential equation with NDSolveValue. The following cell has to be activated first (constants definition, seed definition, mesh definition, ODE definition etc):


Needs["NDSolve`FEM`"];
Needs["NumericalCalculus`"]

(* First define conversion factors between various units*)
CPrec        = 500; (* numerical precission for later calculations*)
constants = 
  SetPrecision[{hbar -> 1.05457*^-34, c0 -> 2.997925*^8, 
    ec -> 1.602177*^-19}, 5 CPrec];
(* Note that ec is the unitless(!!) value of the electric charge, \
which is used to convert eV to J. ev=ec*J *)

RES = SetPrecision[149597870700*m2invMeV, 
   5 CPrec]; (*Earth-Sun distance also known as AE or AU, \
Subscript[r, AU] in the paper*)

aG = SetPrecision[1.30199*10^(-32), 
   5 CPrec]; (*acceleration of earth towards sun in MeV, Pitschmann \
value*)

GN = SetPrecision[6.70861*^-45, 
   5 CPrec]; (*newtons constant in MeV^(-2), Pitschmann value*)
mpl = SetPrecision[2.4353635268*^21, 5 CPrec];(*planck mass in MeV*)
VC = SetPrecision[6.73304*^-34, 
   5 CPrec]; (*in MeV^4, this is the effective potential for the \
vacuum density and the vaccuum expactation value of \[Phi]*)

m2invMeV  = 
  1/(hbar c0/(1*^6*ec)) /. 
   constants; (*to be read as convert meter to MeV^-1*)
kg2MeV       = 1/(ec 10^6/(c0^2)) /. constants; 
s2invMeV  = 1/(hbar/(1*^6 ec)) /. constants; 
N2invMeV2 = kg2MeV m2invMeV/s2invMeV^2; (*Newton to inverse MeV^2*)
kgm32MeV4 = 
  kg2MeV/m2invMeV^3; (*to be read as kg/m^3 converted to MeV^4 / \
density to natural units*)

(*Next define various constants / parameters that are use in the LLR \
analysis*)

rhoV = SetPrecision[1*^-15, 5 CPrec]; (* vacuum density in MeV^4 *)
rhoE = SetPrecision[0.0000236851, 
   5 CPrec]; (*mean density of the earth in MeV^4 *)
rhoS = SetPrecision[1408*kgm32MeV4, 
   5 CPrec]; (* mean density of the sun in MeV^4 *)
rhoM = SetPrecision[0.0000143877, 
   5 CPrec]; (*mean density of the moon in MeV^4, Pitschmannvalue*)
rhoN = SetPrecision[1.435550104*10^10, 
   5 CPrec]; (*neutron density in MeV^4, Pitschmannvalue*)

RE = SetPrecision[6378100*m2invMeV, 
   5 CPrec]; (*equatorial radius of the earth in MeV*)
RS = SetPrecision[695700000*m2invMeV, 
   5 CPrec]; (*equatorial radius of the sun in MeV*)
RM = SetPrecision[1738100*m2invMeV, 
   5 CPrec]; (*equatorial radius of the moon in MeV*)
RN = SetPrecision[0.0025, 
   5 CPrec]; (*Neutron radius MeV^-1 - Pitschmann value*)

ME = SetPrecision[5.9724*10^(24)*kg2MeV, 
   5 CPrec]; (*total mass of the earth in MeV*)
MS = SetPrecision[1.9891*10^(30)*kg2MeV, 
   5 CPrec]; (*total mass of the sun in MeV*)
MN = SetPrecision[939.565346, 
   5 CPrec]; (* neutron mass in MeV, Pitschmann value *)

REM = SetPrecision[385000558.4*m2invMeV, 
   5 CPrec]; (*Earth-Moon distance in MeV^-1*)   
REMax = SetPrecision[405000000*m2invMeV, 
   5 CPrec]; (*Max. Earth-Moon distance in MeV^-1*) 

\[Rho]V = rhoV;
\[Rho]E = rhoE;
\[Rho]S = rhoS;
\[Rho]M = rhoM;
\[Rho]N = rhoN;

\[Rho]M2 = 
  SetPrecision[2514.0`*kgm32MeV4, 5 CPrec]; (*mirror density cannex*)

ddr = SetPrecision[5.0677398985414185`*^7, 1000]; 

phirho[A_, \[Lambda]_, \[Beta]_, \[Rho]_] := 
 SetPrecision[
  mpl/\[Lambda] If[
    2*Log10[\[Lambda]] + \[Beta] - Log10[A] - Log10[\[Rho]] < 1000000,
     ProductLog[(\[Lambda]^2 10^\[Beta])/(A \[Rho])], 
    Log[\[Lambda]^2/(A \[Rho])] + 
     1/(6 (\[Beta] Log[10] + 
            Log[\[Lambda]^2/(A \[Rho])])^3) (6 \[Beta] Log[
          10] (\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])])^3 - 
        6 (-1 + \[Beta] Log[10] + 
           Log[\[Lambda]^2/(A \[Rho])]) (1 + \[Beta]^2 Log[10]^2 + 
           2 \[Beta] Log[10] Log[\[Lambda]^2/(A \[Rho])] + 
           Log[\[Lambda]^2/(A \[Rho])]^2) Log[\[Beta] Log[10] + 
           Log[\[Lambda]^2/(A \[Rho])]] + 
        3 (-3 + \[Beta] Log[10] + 
           Log[\[Lambda]^2/(A \[Rho])]) Log[\[Beta] Log[10] + 
            Log[\[Lambda]^2/(A \[Rho])]]^2 + 
        2 Log[\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])]]^3)], 
  1000] 

murho[\[Lambda]_, \[Rho]_, A2_] := 
  SetPrecision[
   1/mpl Sqrt[\[Lambda]^2 V0[A2, \[Lambda]] Exp[
        SetPrecision[-\[Lambda] phirho[\[Lambda], \[Rho], A2]/mpl, 
         5 CPrec]] + A2 \[Rho]], 5 CPrec];

Veff[\[Lambda]_, \[Rho]_, \[Phi]_, A2_] := 
  SetPrecision[
   V0[A2, \[Lambda]] Exp[
      SetPrecision[-\[Lambda] \[Phi]/mpl, 
       5 CPrec]] + (A2 \[Rho])/(2 mpl^2) \[Phi]^2, 5 CPrec];

j[x_] := 3/x^2 (1 - Tanh[x]/x)

js[x_] := 
 1 - (2 x^2)/5 + (17 x^4)/105 - (62 x^6)/945 + (1382 x^8)/51975

SC[\[Lambda]_, \[Rho]_, A2_, R_] := 
 If[murho[\[Lambda], \[Rho], A2] R <= 0.01, 
  SetPrecision[
   js[murho[\[Lambda], \[Rho], A2] R]/(1 + 
      murho[\[Lambda], rhoV, A2]/murho[\[Lambda], \[Rho], A2] Tanh[
        murho[\[Lambda], \[Rho], A2] R]), 5 CPrec], 
  SetPrecision[
   j[murho[\[Lambda], \[Rho], A2] R]/(1 + 
      murho[\[Lambda], rhoV, A2]/murho[\[Lambda], \[Rho], A2] Tanh[
        murho[\[Lambda], \[Rho], A2] R]), 
   5 CPrec]] (*screening charge*)

Q[\[Lambda]_, \[Rho]_, A2_, R_] := SC[\[Lambda], \[Rho], A2, R]

\[Mu][A_, \[Lambda]_, \[Beta]_, \[Rho]_] := 
 SetPrecision[ 
  Sqrt[SetPrecision[\[Lambda]^2  E^
        SetPrecision[-\[Lambda] \[Phi]M[
             A, \[Lambda], \[Beta], \[Rho]]/mpl + \[Beta]/Log10[E], 
         1000] + A \[Rho], 1000]]/mpl, 1000]

\[Phi]M[A_, \[Lambda]_, \[Beta]_, \[Rho]_] := 
 phirho[A, \[Lambda], \[Beta], \[Rho]]

\[Mu]0[\[Lambda]_?NumberQ, 
  A2_?NumberQ, \[Beta]_?NumberQ, \[Phi]0_?NumberQ, \[Rho]V2_?
   NumberQ] := 
 SetPrecision[
  1/mpl Sqrt[\[Lambda]^2 E^
       SetPrecision[-\[Lambda] \[Phi]0/mpl + \[Beta]/Log10[E], 1000] +
      A2 \[Rho]V2], 1000]

D0[\[Lambda]_?NumberQ, 
  A2_?NumberQ, \[Beta]_?NumberQ, \[Phi]0_?NumberQ, \[Rho]V2_?
   NumberQ] := 
 SetPrecision[\[Lambda] /mpl E^
     SetPrecision[-\[Lambda] \[Phi]0/mpl + \[Beta]/Log10[E], 
      1000] - (A2 \[Rho]V2)/mpl^2 \[Phi]0, 1000]

\[Phi]d[\[Lambda]_, A2_, \[Beta]_, \[Phi]0_, d_, \[Rho]V2_] := 
 SetPrecision[\[Phi]0 + 
   D0[\[Lambda], 
      A2, \[Beta], \[Phi]0, \[Rho]V2]/\[Mu]0[\[Lambda], 
       A2, \[Beta], \[Phi]0, \[Rho]V2]^2 (1 - 
      Cosh[SetPrecision[\[Mu]0[\[Lambda], 
          A2, \[Beta], \[Phi]0, \[Rho]V2]*d, 5 CPrec]]), 5 CPrec]

\[Phi]2[\[Lambda]_, A2_, \[Beta]_, \[Phi]0_, d_, 
  z_, \[Rho]M2_, \[Rho]V2_] := 
 SetPrecision[
  Piecewise[{{\[Phi]0 + 
      D0[\[Lambda], 
         A2, \[Beta], \[Phi]0, \[Rho]V2]/\[Mu]0[\[Lambda], 
          A2, \[Beta], \[Phi]0, \[Rho]V2]^2 (1 - 
         Cosh[SetPrecision[\[Mu]0[\[Lambda], 
             A2, \[Beta], \[Phi]0, \[Rho]V2]*(z), 5 CPrec]]), 
     d - Abs[z] > 0},
    {\[Phi]M[
       A2, \[Lambda], \[Beta], \[Rho]M2] + (\[Phi]d[\[Lambda], 
          A2, \[Beta], \[Phi]0, d, \[Rho]V2] - \[Phi]M[
          A2, \[Lambda], \[Beta], \[Rho]M2]) E^
        SetPrecision[-\[Mu][
            A2, \[Lambda], \[Beta], \[Rho]M2] (Abs[z] - d), 5 CPrec], 
     Abs[z] - d > 0}}], 5 CPrec]

\[Phi]\[Rho][A_, \[Lambda]_, \[Beta]_, \[Rho]_] := 
 SetPrecision[
  mpl/\[Lambda] If[
    2*Log10[\[Lambda]] + \[Beta] - Log10[A] - Log10[\[Rho]] < 1000000,
     ProductLog[(\[Lambda]^2 10^\[Beta])/(A \[Rho])], 
    Log[\[Lambda]^2/(A \[Rho])] + 
     1/(6 (\[Beta] Log[10] + 
            Log[\[Lambda]^2/(A \[Rho])])^3) (6 \[Beta] Log[
          10] (\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])])^3 - 
        6 (-1 + \[Beta] Log[10] + 
           Log[\[Lambda]^2/(A \[Rho])]) (1 + \[Beta]^2 Log[10]^2 + 
           2 \[Beta] Log[10] Log[\[Lambda]^2/(A \[Rho])] + 
           Log[\[Lambda]^2/(A \[Rho])]^2) Log[\[Beta] Log[10] + 
           Log[\[Lambda]^2/(A \[Rho])]] + 
        3 (-3 + \[Beta] Log[10] + 
           Log[\[Lambda]^2/(A \[Rho])]) Log[\[Beta] Log[10] + 
            Log[\[Lambda]^2/(A \[Rho])]]^2 + 
        2 Log[\[Beta] Log[10] + Log[\[Lambda]^2/(A \[Rho])]]^3)], 1000]
EQN\[Phi]0[\[Lambda]_, A2_, \[Beta]_, \[Phi]0_, \[Rho]_, \[Rho]V2_, 
  d_] := SetPrecision[-\[Phi]0 + ((mpl \[Lambda] - 
        A2 E^SetPrecision[(\[Lambda] \[Phi]0)/mpl - \[Beta] Log[10], 
           1000] \[Rho]V2 \[Phi]0) (-1 + 
        Cosh[SetPrecision[(d Sqrt[
              E^SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[
                    10], 1000] \[Lambda]^2 + A2 \[Rho]V2])/mpl, 
          1000]]))/(\[Lambda]^2 + 
      A2 E^SetPrecision[(\[Lambda] \[Phi]0)/mpl - \[Beta] Log[10], 
         1000] \[Rho]V2) + \[Phi]\[Rho][
    A2, \[Lambda], \[Beta] , \[Rho]] + ((E^
          SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[10], 
           1000] mpl \[Lambda] - 
        A2 \[Rho]V2 \[Phi]0) Sinh[(d Sqrt[
           E^SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[
                  10], 1000] \[Lambda]^2 + A2 \[Rho]V2])/mpl])/(Sqrt[
       E^SetPrecision[-((\[Lambda] \[Phi]0)/mpl) + \[Beta] Log[10], 
           1000] \[Lambda]^2 + A2 \[Rho]V2] Sqrt[
       A2 \[Rho] (1 + (\[Lambda] \[Phi]\[Rho][
              A2, \[Lambda], \[Beta] , \[Rho]])/mpl)]), 1000]
\[Phi]0[\[Lambda]_?NumericQ, A2_, \[Beta]_, d_, \[Rho]_, \[Rho]V2_] :=
   Block[{res, startp},
   startp = 
    Abs[SetPrecision[\[Phi]\[Rho][A2, \[Lambda], \[Beta], \[Rho]V2], 
      5 CPrec]];
   res = 
    N[(phi0 /. 
       SetPrecision[
        FindRoot[
         Re[SetPrecision[
           EQN\[Phi]0[SetPrecision[\[Lambda], 1000], 
            SetPrecision[A2, 1000], SetPrecision[\[Beta], 1000], 
            SetPrecision[phi0, 1000], SetPrecision[\[Rho], 1000], 
            SetPrecision[\[Rho]V2, 1000], SetPrecision[d, 1000]], 
           5 CPrec]], {phi0, startp, 0, 
          SetPrecision[1.5, 5 CPrec] startp}, AccuracyGoal -> 5 CPrec,
          PrecisionGoal -> 5 CPrec, WorkingPrecision -> 5 CPrec], 
        5 CPrec]), 5 CPrec];
   test = checkval[res];
   If[test != NaN,
    If[Im[test]/Re[test] > 1*^-5, 
     Print["Warning: getphi0i resulted in non-negligible Complex \
value at \[Lambda]=", N[\[Lambda], 6], " and A2=", N[A2, 6]]]];
   Abs[res]];(*10 \[Mu]m in MeV*)

(*Run entire cell*)
Clear[ResIncrease, s, r, z, domb];

Clear[\[Rho]V, \[Rho]M]

\[Rho]V := 
  2.282744303895228569942873684140791998747198853011172394898362081328\
00265957484953105449676513671875`1000.*^-20;

rhoV = \[Rho]V;

\[Rho]M = SetPrecision[2514.0`*kgm32MeV4, 5 CPrec];

rhoM = \[Rho]M;

ResIncrease = 10^0;
s = 1/ResIncrease;

r[z_] :=  m Tan[s*z];
z[r_] := ArcTan[r/m]/s;

m = m2invMeV;
dd = SetPrecision[ArcTan[10.0`*^-6]/s, 5 CPrec]; (* plate separation *)
dp = SetPrecision[ArcTan[100.0`*^-6]/s, 
   5 CPrec]; (*thickness of upper plate*)
ddr = 10.0`*^-6*m2invMeV;
dpr = 100.0`*^-6*m2invMeV


initialval1[r_, 
  A2_?NumberQ, \[Lambda]_?NumberQ, \[Beta]_?
   NumberQ] := \[Phi]2[\[Lambda], 
  A2, \[Beta], \[Phi]0[\[Lambda], A2, \[Beta], 
   ddr/2, \[Rho]M, \[Rho]V], ddr/2, r, \[Rho]M, \[Rho]V]


Seed[z_, A2_?NumberQ, \[Lambda]_?NumberQ, \[Beta]_?NumberQ] := 
  initialval1[r[z], A2, \[Lambda], \[Beta]];


Clear[rho];
rho[z_] := 
  SetPrecision[
   Piecewise[{{\[Rho]V, (-0.5*dd < z < 0.5*dd)}, {\[Rho]M, 
      z <= -0.5*dd || z >= -0.5*dd}}], 5 CPrec];
\[Rho][z_] := rho[z]

Clear[C1, C2, C3, C4];

C1[z_?NumberQ] := SetPrecision[mpl/(m^2*s^2) Cos[s*z]^4, 5 CPrec];

C2[z_?NumberQ] := 
  SetPrecision[-2 mpl/(m^2*s) Cos[s*z]^4 Tan[s*z], 5 CPrec];

C3[z_?NumberQ, A2_?NumberQ] := 
  SetPrecision[-A2/mpl \[Rho][z], 5 CPrec];

C4[f_, \[Lambda]_?NumberQ, A2_?NumberQ, \[Beta]_?NumberQ] := 
 If[- ((\[Lambda] f)/mpl) + Log[\[Lambda]] + Log[10]*\[Beta] < 0, 
  SetPrecision[E^
   SetPrecision[-Sqrt[(- ((\[Lambda] f)/mpl) + Log[\[Lambda]] + 
       Log[10]*\[Beta])^2], 5 CPrec], 5 CPrec], 
  SetPrecision[E^
   SetPrecision[
    Sqrt[(- ((\[Lambda] f)/mpl) + Log[\[Lambda]] + 
      Log[10]*\[Beta])^2], 5 CPrec], 5 CPrec]]


Residual[z_, A2_, \[Lambda]_, \[Beta]_, u_] := 
 Abs[C1[z]*u''[z] + C2[z]*u'[z] + C3[z, A2]*u[z] + 
   C4[u[z], \[Lambda], A2, \[Beta]]]
RelResidual[z_, A2_, \[Lambda]_, \[Beta]_, u_] := 
 Residual[z, A2, \[Lambda], \[Beta], u]/
  Max[Abs[C1[z]*u''[z]], Abs[C2[z]*u'[z]], Abs[C3[z, A2]*u[z]], 
   Abs[C4[u[z], \[Lambda], A2, \[Beta]]]]

coordlist[dmin_, ddelta_, dfactor_, nmin_, nmax_] := 
  Block[{i, dummy},
   dummy = Table[N[dmin, 5 CPrec], nmax - nmin + 1];
   If[nmin == 0,
    For[i = 1, i <= nmax, i++,
      dummy[[i + 1]] = dummy[[i]] + N[ddelta*dfactor^i, 5 CPrec];
      ];,
    dummy[[1]] += Sum[N[ddelta*dfactor^i, 5 CPrec], {i, 1, nmin}];
    For[i = 2, i <= nmax - nmin + 1, i++,
     dummy[[i]] = 
       dummy[[i - 1]] + N[ddelta*dfactor^(nmin + i - 1), 5 CPrec];
     ];];
   N[dummy, 5 CPrec]
   ];


Clear[bvp, x, \[Lambda], A2, \[Beta]];
bvp = {SetPrecision[
      C1[z]*D[D[u[z], z], z] + C2[z]*D[u[z], z] + C3[z, A2]*u[z] + 
       C4[u[z], \[Lambda], A2, \[Beta]], 5 CPrec] == 
     SetPrecision[0, 5 CPrec],
    SetPrecision[u[SetPrecision[-\[Pi]/2, 5 CPrec]], 5 CPrec] == 
     N[SetPrecision[\[Phi]M[A2, \[Lambda], \[Beta], \[Rho]M], 
       5 CPrec], 5 CPrec],
    SetPrecision[u[SetPrecision[\[Pi]/(2*s), 5 CPrec]], 5 CPrec] == 
     SetPrecision[
      N[\[Phi]M[A2, \[Lambda], \[Beta], \[Rho]M], 5 CPrec], 
      5 CPrec]} /. m -> m2invMeV;
dist1 = SetPrecision[\[Pi]/(2.0`*s), 
   5 CPrec];(* from -\[Infinity] to the border of the lower plate*)


(* Was ist dist 2? dp/2 nicht legitim?*)

dist2 = SetPrecision[dp/2., 5 CPrec];(* half thickness of the plate*)

(*was sollen die kommenden Größen sein?*)
nMaxRes = 110*^5;(*max. resolution in the gap 110*^6*)
npts1 = 80;(*points in free vacuum and free material, 400*)
npts2 = 80;(*points in the half-gap, 500*)
npts3 = 80;(*points in the upper plate, 500*)
nFinePoints = 30;(*120*)
minMeshSize = SetPrecision[(s*dd/nMaxRes), 5 CPrec];

f1 =(*von -\[Pi]/2 bis -dd/2*)
  SetPrecision[(\[Delta] /. 
     Flatten[FindRoot[
        SetPrecision[
          Sum[N[minMeshSize, 5 CPrec]*\[Delta]^N[n, 5 CPrec], {n, 1, 
            npts1 - 1}], 5 CPrec] == 
         SetPrecision[
          dist1 - dd/2 - N[nFinePoints*minMeshSize, 5 CPrec], 
          5 CPrec], {\[Delta], 1.01`}, MaxIterations -> 3000, 
        PrecisionGoal -> 5 CPrec + 2, AccuracyGoal -> 5 CPrec + 2, 
        WorkingPrecision -> 5 CPrec]][[1]]), 5 CPrec];
(*f2= (*half thickness of the upper \
plate*)SetPrecision[(\[Delta]/.Flatten[FindRoot[SetPrecision[Sum[N[\
minMeshSize,5 CPrec]*\[Delta]^N[n,5 CPrec],{n,1,npts3-1}],5 CPrec]\
\[Equal]SetPrecision[dist2-N[nFinePoints*minMeshSize,5 CPrec],5 \
CPrec],{\[Delta],1.01`},MaxIterations\[Rule]3000,PrecisionGoal\[Rule]\
5 CPrec+2,AccuracyGoal\[Rule]5 CPrec+2,WorkingPrecision\[Rule]5 \
CPrec]][[1]]),5CPrec];*)
(*f3=(*points above upper plate divided by \
2*)SetPrecision[(\[Delta]/.Flatten[FindRoot[SetPrecision[Sum[N[\
minMeshSize,5 CPrec]*\[Delta]^N[n,5 CPrec],{n,1,npts1-1}],5 CPrec]\
\[Equal]SetPrecision[0.5*(dist1-dp-dd/2)-N[nFinePoints*minMeshSize,5 \
CPrec],5 \
CPrec],{\[Delta],1.01`},MaxIterations\[Rule]3000,PrecisionGoal\[Rule]\
5 CPrec+2,AccuracyGoal\[Rule]5 CPrec+2,WorkingPrecision\[Rule]5 \
CPrec]][[1]]),5CPrec];*)
f4 =(*half gab between plates*)
  SetPrecision[(\[Delta] /. 
     Flatten[FindRoot[
        SetPrecision[
          Sum[N[minMeshSize, 5 CPrec]*\[Delta]^N[n, 5 CPrec], {n, 1, 
            npts2 - 1}], 5 CPrec] == 
         SetPrecision[dd/2 - N[nFinePoints*minMeshSize, 5 CPrec], 
          5 CPrec], {\[Delta], 1.001`}, MaxIterations -> 3000, 
        PrecisionGoal -> 5 CPrec + 2, AccuracyGoal -> 5 CPrec + 2, 
        WorkingPrecision -> 5 CPrec]][[1]]), 5 CPrec];

points = SetPrecision[Sort[N[Join[
      (* lower plate *)
      coordlist[
       N[-nFinePoints*minMeshSize - dd/2, 5 CPrec], -minMeshSize, f1, 
       1, npts1 - 1],
      Table[
       N[-i*minMeshSize - dd/2, 5 CPrec], {i, 0, nFinePoints - 1}],
      Table[N[i*minMeshSize - dd/2, 5 CPrec], {i, 1, nFinePoints}],
      Table[-\[Pi]/2 + N[i*minMeshSize, 5 CPrec], {i, 1, 
        nFinePoints}] ,
      
      (*in between plates*)
      coordlist[N[nFinePoints*minMeshSize - dd/2, 5 CPrec], 
       minMeshSize, f4, 1, npts2 - 2], 
      coordlist[
       dd/2 - N[nFinePoints*minMeshSize, 5 CPrec], -minMeshSize, f4, 
       1, npts2 - 1],
      Table[
       dd/2 - N[i*minMeshSize, 5 CPrec], {i, 0, nFinePoints - 1}](* 
      inter-spacing *),
      Table[dd/2 + N[i*minMeshSize, 5 CPrec], {i, 1, nFinePoints - 1}],
      
      (* upper plate *)
      coordlist[N[nFinePoints*minMeshSize + dd/2, 5 CPrec], 
       minMeshSize, f1, 1, npts1 - 1],
      Table[\[Pi]/2 - N[i*minMeshSize, 5 CPrec], {i, 1, nFinePoints}]
      
      ]]](*free space from right*), 5 CPrec];
points[[1]] = -\[Pi]/2;
points[[-1]] = \[Pi]/(2*s);

domb = SetPrecision[
   ToElementMesh["Coordinates" -> Transpose[{points}], 
    "MeshElements" -> {LineElement[
       Table[{i, i + 1}, {i, 1, Length[points] - 1}]]}, 
    "BoundaryElements" -> {(*external boundaries:*)
      PointElement[{{1}, {Length[points]}}]}], 5 CPrec];
Length[points]        

After the above code has been activated one can simply solve the ODE by deciding on the parameters (following code). This gives the correct result for [Lambda]=10^(27), not for [Lambda]=10^(31)

Clear[\[Lambda]t, A2, \[Lambda], A2, Y1, Y2, \[Beta]]
\[Lambda] = 10^31;
A2 = 10^40;
\[Beta] = 1;
Y1 = NDSolveValue[Evaluate[bvp], u, z \[Element] domb, 
  InitialSeeding -> {u[z] == Seed[z, A2, \[Lambda], \[Beta]]}]
Plot[{Y1[z], Seed[z, A2, \[Lambda], \[Beta]]}, {z, -dd, dd}]
    

Any help or insight would be appreciated, thank you!

$\endgroup$
1
  • 1
    $\begingroup$ The finite element method is only available for machine precision. So setting precision for FEM does not do anything. Next, perhaps the ODE is strongly nonlinear, if that is the case you can use ParametricNDSolve for a solver restart. The SolidMechanics monograph has an example of how to deal with that here. $\endgroup$
    – user21
    Commented Jan 10, 2023 at 12:40

1 Answer 1

1
$\begingroup$

I did find a solution to the above problem, and I am going to post an answer in case somebody comes across similiar convergence issues. For some parameters (e.g. [beta]=1, A2=10^(40), \lambda=10^(31) ) mathematica does return a solution after only few iterations (without error messages), but the solution is not converged.

If I use the returned solution as a seed to a new simulation, Mathematica returns a slightly different solution after only few iterations. Repeating this procedure several times eventually returns the correct and fully converged solution. I would like to note that for some reason mathematica never goes through the full 100 iterations, even though it could converge to the correct solution with only a few more iterations. I suppose that his is a bug..

$\endgroup$

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.