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$ Jan 1 at 12:37
  • $\begingroup$ Well, l2num should be very small, for instance, of 10^-16 order. $\endgroup$ Jan 1 at 13:16
4
$\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$ Jan 2 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$ Jan 2 at 12:16
  • $\begingroup$ In fact, NIntegrate[Evaluate[s1[eps]], {eps, eps0, l1numFR}] == 4688.49. What am I doing wrong? $\endgroup$ Jan 2 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
    Jan 2 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
    Jan 4 at 0:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.