0

I am trying to solve for the limits of integration using a piecewise function. I'm doing this in a loop over a list of data, mostly because I'm not well-versed enough in Mathematica to make this a snazzy function or use Map. I've had some success doing the following with an analytical function (i.e. without a piecewise function) and now I'm getting some warnings and it's taking much longer than expected:

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

Here is the code I am running:

(* Create some data *)
radTh = RandomReal[{100, 2000}, {10, 3}];
radTh[[;; , 3]] = {320.128`, 328.839`, 317.76`, 321.142`, 323.156`, 
                  325.015`, 323.975`, 315.921`, 318.175`, 315.006`};

rhoS = radTh;
rhoS[[;; , 3]] = 158.551;

radBd = radTh;
radBd[[;; , 3]] = {-259.5056567581`, -270.5568561712`, -261.82074464040005`, 
  -267.54610008360004`, -271.9701557155`, -276.1683057652`, -276.603896472`, 
  -266.94477017400004`, -267.56731240690004`, -264.792137731`};

rT = radTh;
rT[[;; , 3]] = {3.7904278208442067`*^-6, 3.893569116661423`*^-6, 3.762389870150237`*^-6, 
  3.802433936555222`*^-6, 3.826280403066056`*^-6, 3.8482916152029184`*^-6, 
  3.835977650371107`*^-6, 3.740615465029371`*^-6,    3.7673036157321617`*^-6, 
  3.7297815440475374`*^-6};

rho[rhos_, Z_, Zp_] := 917. - (917. - rhos)*Exp[(-Z/Zp)]

radSf = radTh;
radSf[[;; , 3]] = {60.62234324190001`, 58.282143828799974`, 55.93925535959994`, 
  53.59589991639996`, 51.18584428449998`, 48.8466942348`, 47.37110352800005`, 
  48.97622982599995`, 50.60768759309997`, 50.21386226899995`};

c2 = 845.*10^(-6);
c = 299792458.;
g = 9.8;

RSRthick = radTh;
RSRwt = radTh;
HSrsr = radTh;
For[i = 1, i <= 10, i++,

  tempLinearProfile = {{0., rhoS[[i, 3]]}, {9.9, 
  rho[485.6, 9.9, 48.]}};
  tempLinearProfileInterp = Interpolation[tempLinearProfile, InterpolationOrder -> 1];

  LDF[Z_] = Piecewise[{{tempLinearProfileInterp[Z], Z <= 9.9}, 
                       {rho[485.6, Z, 48.], Z > 9.9}}];

  RSRthick[[i, 3]] =  Hrsr /. First@
   FindRoot[rT[[i, 3]] == (2. *Integrate[(LDF[Z]*c2 + 1.), {Z, 0., Hrsr}, 
     Assumptions -> Hrsr > 0])/c , {Hrsr, radTh[[i, 3]]}];

  RSRwt[[i, 3]] = NIntegrate[LDF[Z]*g, {Z, 0., RSRthick[[i, 3]]}];

  HSrsr[[i, 3]] = -Drsr /. First@FindRoot[ 
    Drsr*1028. - 
    Integrate[LDF[Z], {Z, radSf[[i, 3]], (radSf[[i, 3]] + Drsr)},
         Assumptions -> Drsr > 0] == 
    Integrate[LDF[Z], {Z, 0., radSf[[i, 3]]}, Assumptions -> Drsr > 0], 
    {Drsr, -radBd[[i, 3]]}]

] // Timing

When I used an analytic function for LDZ I had some luck using NSolve (with rationalize) but along the way switched to FindRoot. I'm not sure what the differences are between these methods? What is more appropriate for solving for the limit of integration for analytic function vs piecewise functions? I imagine that having a continuous derivative plays into the different methods that are used in Mathematica, but I'm just guessing.

It would be great to gain some input on the difference between NSolve and FindRoot methods, best practices for implementing this type of routine/solves in Mathematica, appropriate "line search" step sizes, and ways to speed up the calculations (my data is ~100k points). Thanks!

user38307
  • 41
  • 3
  • 2
    A MINIMAL example would be much better – Dr. belisarius Mar 07 '16 at 18:13
  • Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the [faq]! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Dr. belisarius Mar 07 '16 at 18:13
  • I'm voting to close this question as off-topic because it is too localized; i.e, it applies only to the local situation and needs of its poster and answers will not benefit others. – m_goldberg Mar 17 '16 at 14:37
  • Ok, no problem. I think my issue could be of broader appeal if reasked in a way that was pointing to the issues that are sort of addressed in http://mathematica.stackexchange.com/questions/7924/alternatives-to-procedural-loops-and-iterating-over-lists-in-mathematica – user38307 Mar 18 '16 at 22:06
  • Ok, no problem @m_goldberg,@Dr. belisarius. The way I'm doing this sort of works.

    I'm a bit lost on avoiding using For loops, but I can ask that as an independent post that may have broader appeal. I've looked through the "common pitfalls" and "alternatives to procedural loops" links.

    I'm assuming that avoiding a for loop has a performance benefit. What I'm unclear on is what is the appropriate approach for running multiple statements that are dependent on each other (i.e. can't be run independently).

     – user38307 Mar 18 '16 at 22:27

0 Answers0