Find right end-point satisfying an integral constraint

I am solving the following system of first-order ODEs with a variable right-end point l1num and boundary condition at that point n[l1num] == l2num:

a1 = 5.85476*10^8;
a2 = 1.24911*10^45;
a3 = 9.648*10^46;

eps0 = 10^-5;

epsilon = {0, 0.03, 0.05, 0.08, 0.10, 0.13, 0.15, 0.18, 0.20, 0.23,
0.25, 0.28, 0.30, 0.38, 0.45, 0.55, 0.68, 0.83, 1.00, 1.98, 2.00,
2.10, 4.23, 5.10, 5.95, 6.00, 6.05, 7.95, 8.00, 8.05, 10.10, 12.20,
15.20, 18.20, 23.10, 26.20, 30.90, 37.20, 45.60, 51.00, 63.00,
79.00, 80.00, 81.00, 94.00, 95.00, 96.00, 118.00, 149.00, 150.00,
151.00, 188.00, 231.00, 284.00, 333.00, 396.00, 483, 584, 715, 895,
1000};
sigma = {1.26 10^-19, 2.82 10^-20, 1.41 10^-20, 8.71 10^-21,
6.21 10^-21, 4.65 10^-21, 3.68 10^-21, 3.02 10^-21, 2.54 10^-21,
2.19 10^-21, 1.91 10^-21, 1.69 10^-21, 1.51 10^-21, 2.27 10^-21,
3.17 10^-21, 4.58 10^-21, 6.67 10^-21, 9.63 10^-21, 1.37 10^-20,
2.67 10^-20, 2.74 10^-20, 2.85 10^-20, 6.10 10^-20, 7.65 10^-20,
8.93 10^-20, 9.30 10^-20, 9.08 10^-20, 1.19 10^-19, 1.22 10^-19,
1.21 10^-19, 1.40 10^-19, 1.38 10^-19, 1.23 10^-19, 1.02 10^-19,
7.66 10^-20, 6.54 10^-20, 5.42 10^-20, 4.43 10^-20, 3.51 10^-20,
3.09 10^-20, 2.53 10^-20, 2.02 10^-20, 2.19 10^-20, 1.98 10^-20,
1.72 10^-20, 1.74 10^-20, 1.68 10^-20, 1.38 10^-20, 1.11 10^-20,
1.14 10^-20, 1.09 10^-20, 9.15 10^-21, 7.63 10^-21, 6.32 10^-21,
5.46 10^-21, 4.64 E 10^-21, 3.83 10^-21, 3.18 10^-21, 2.55 10^-21,
1.95 10^-21, 1.70 10^-21};
interpolant = Transpose@Join[{epsilon}, {sigma}];
sigma1 = Interpolation[interpolant, InterpolationOrder -> 1];

solODE = NDSolve[{
n1'[eps] == n[eps],
a1 (n[eps] - 2 eps n'[eps]) + eps sigma1[eps]^2 (a2 (n[eps] - 2 eps n'[eps]) - a3 eps n[eps]) == 0,
n1[l1num] == 1, n[l1num] == l2num}, {n, n1}, {eps, eps0, l1num}];


The solution of which must satisfy the integral constraint

NIntegrate[n[eps], {eps, eps0, l1num}] == 1.


I tried to fix l2num and store the solution in a table with respect to l1num:

l2num = 10^-16;
solODE2[l1_] := NDSolve[{
n1'[eps] == n[eps],
a1 (n[eps] - 2 eps n'[eps]) + eps sigma1[eps]^2 (a2 (n[eps] - 2 eps n'[eps]) - a3 eps n[eps]) == 0,
n1[l1] == 1, n[l1] == l2num}, {n, n1}, {eps, eps0, l1}];

tabl1l2 = Table[Abs[1 - NIntegrate[Evaluate[n[eps] /. solODE2[l1]], {eps, eps0, l1}]], {l1, 10, 20, .001}];


Then,

Min[tabl1l2] = 0.745373


Thus, in this way, I can get some very rough solution, but it takes a long time since I need a finer division in table.

I also tried to locate the required value using FindRoot':

FindRoot[1 - NIntegrate[Evaluate[n[eps] /. solODE2[l1]], {eps, eps0, l1}], {l1, 10}]


which results in

NDSolve::ndsv: Cannot find starting value for the variable n.


Any idea is highly appreciated.

Update

I think, it will be more appropriate to change the problem statement as follows. How to solve the following ODE

solODEred = NDSolve[{a1 (n[eps] - 2 eps n'[eps]) + eps sigma1[eps]^2 (a2 (n[eps] - 2 eps n'[eps]) - a3 eps n[eps]) == 0,
n[l1num] == l2num}, {n, n1}, {eps, eps0, l1num}];


such that its solution satisfies the integral constraint:

NIntegrate[n[eps], {eps, eps0, l1num}] == 1.


We can fix l2num = 10^-16. We can even fix l1num (that would be even better). Apparently, the addition of a new equation n1[eps] == n[eps] and corresponding condition n1[l1num] == 1 is not of much help.

• Differentiating the integral constraint gives D[Integrate[n[eps], {eps, eps0, l1num}], l1num]==n[l1num]==0 , perhaps in contradiction to n[l1num]==l2num? Commented Jan 1, 2022 at 12:37
• Well, l2num should be very small, for instance, of 10^-16 order. Commented Jan 1, 2022 at 13:16

With the constants defined as in the question, compute

{s1, s2} = NDSolveValue[{n1'[eps] == n[eps], a1 (n[eps] - 2 eps n'[eps]) +
eps sigma1[eps]^2 (a2 (n[eps] - 2 eps n'[eps]) - a3 eps n[eps]) == 0,
n1[eps0] == 0, n[eps0] == 1}, {n, n1}, {eps, eps0, 12}];


which differs from the computation given in the question by the boundary conditions, n1[eps0] == 0, n[eps0] == 1 and the upper bound on the integration. Because the ODEs are linear, this computation determines the desired values of n and n1 (which is the integral of n)) up to the same multiplicative constant. In other words, the computation above gives the desired value for s1[eps]/s2[eps], and s1[l1num]/s2[l1num] == l2num determines l1num. Ploting

LogPlot[s1[eps]/s2[eps], {eps, eps0, 12}, ImageSize -> Large,
AxesLabel -> {eps, n/n1}, LabelStyle -> {15, Bold, Black}]


we find that the result has inadequate precision. To eliminate this difficulty, Rationalize all constants (e.g., a1 = Rationalize[5.85476*10^8, 0];, include the option, WorkingPrecision -> 30 in NDSolve, and repeat the computation to obtain

With the added precision, we obtain the desired

FindRoot[s1[l1num]/s2[l1num] == l2num, {l1num, 8}]
(* {l1num -> 10.0917} *)


Addendum: More direct but much slower solutions

Defining solODE2 as in the question but with WorkingPrecision -> 40, MaxSteps -> 30000 allows more direct but much slower computations.

FindRoot[NIntegrate[n[eps] /. solODE2[l1], {eps, eps0, l1},
WorkingPrecision -> 30, PrecisionGoal -> 6] == 1, {l1,8}, Evaluated -> False]


or, more simply but still slow,

FindRoot[n1[eps0] /. solODE2[l1], {l1, 8}, Evaluated -> False]


both yield the earlier result, {l1num -> 10.0917}.

• I appreciate your reply. So, I run the code l2num = 10^-16; {s1, s2} = NDSolveValue[{n1'[eps] == n[eps], a1 (n[eps] - 2 eps n'[eps]) + eps sigma1[eps]^2 (a2 (n[eps] - 2 eps n'[eps]) - a3 eps n[eps]) == 0, n1[eps0] == 0, n[eps0] == 1}, {n, n1}, {eps, eps0, 12}, WorkingPrecision -> MachinePrecision]; l1numFR = l1num /. FindRoot[s1[l1num]/s2[l1num] == l2num, {l1num, 8}]; and it indeed finds a root. Commented Jan 2, 2022 at 12:13
• However, when I enter the root into the initial system: l2num = 10^-16; solODE2[l1_] := NDSolve[{n1'[eps] == n[eps], a1 (n[eps] - 2 eps n'[eps]) + eps sigma1[eps]^2 (a2 (n[eps] - 2 eps n'[eps]) - a3 eps n[eps]) == 0, n1[l1] == 1, n[l1] == l2num}, {n, n1}, {eps, eps0, l1}, WorkingPrecision -> MachinePrecision]; the solution does not satisfy the integral constraint. Commented Jan 2, 2022 at 12:16
• In fact, NIntegrate[Evaluate[s1[eps]], {eps, eps0, l1numFR}] == 4688.49. What am I doing wrong? Commented Jan 2, 2022 at 12:16
• Don't use WorkingPrecision -> MachinePrecision. The precision is inadequate, as shown in my first plot. Rationalize all constants and use WorkingPrecision -> 30. Commented Jan 2, 2022 at 15:37
• Plotting s1[eps]/s2[eps] shows that the answer must be around 6, so use that as a guess. FindRoot[s1[eps]/s2[eps] == l2num, {eps, 6}] yields {eps -> 5.97056}. Always plot the curve you are working with to see what a good initial guess would be. Also, note that the plot is nonsense for eps > 7 due to insufficient precision. Unfortunately, increasing WorkingPrecision beyond about 45` often leads to internal errors that are beyond the user's control. Commented Jan 4, 2022 at 0:18