1
$\begingroup$

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.

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

1 Answer 1

5
$\begingroup$

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}]

enter image description here

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

enter image description here

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}.

$\endgroup$
13
  • $\begingroup$ 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. $\endgroup$ Commented Jan 2, 2022 at 12:13
  • $\begingroup$ 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. $\endgroup$ Commented Jan 2, 2022 at 12:16
  • $\begingroup$ In fact, NIntegrate[Evaluate[s1[eps]], {eps, eps0, l1numFR}] == 4688.49. What am I doing wrong? $\endgroup$ Commented Jan 2, 2022 at 12:16
  • $\begingroup$ Don't use WorkingPrecision -> MachinePrecision. The precision is inadequate, as shown in my first plot. Rationalize all constants and use WorkingPrecision -> 30. $\endgroup$
    – bbgodfrey
    Commented Jan 2, 2022 at 15:37
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented Jan 4, 2022 at 0:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.