Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am studying a set of functions (many of which I know only as a definite integral) and I have assembled into a list. Here is a sample:

Clear[functionsToPlot, totalFunctions]
functionsToPlot[x_] := {
  (*1*) NIntegrate[0.3/(z + 2) Sqrt[z^2 + 2.0 z + 0.7]/((z - 1)^2 + 1.1), {z, 0, x}],
  (*2*) NIntegrate[.03/(z + 2) Sqrt[z^2 + 2.5 z + 1.50 x]/((z - 1)^2 + .03), {z, 0, x}],
  (*3*) 0.5*Sqrt[x + 0.5]/x   
};

What I need is to make a plot of each curve divided by the total:

totalFunctions[x_] := Total[functionsToPlot[x]]

Plot[Evaluate[functionsToPlot[x]/totalFunctions[x]], {x, 0, 3}]

enter image description here

But the problem is that it takes forever. I'm certain the Plot routine is doing the numerical integrals twice (once in the numerator, and again for the denominator). Is there a way to teach Mathematica to form the total as it computes the various curves in the numerator as it produces the plot in realtime (or to make the task run more efficiently)?

share|improve this question

1 Answer 1

up vote 10 down vote accepted

Update

Almost ten times faster again, or about 90 times faster than the OP's way (0.069 sec v. 5.46 sec):

For the second integral, we can find its derivative with respect to x and then integrate with NDSolve. The derivative of the integral has two components, one from differentiating under the integral sign dxdz1 and one from plugging in the limit of integration. There is a coincidence here, that sometimes happens, namely, that the derivative with respect to x of the integrand can be integrated with respect to z relatively easily.

ad = Integrate[
   D[Rationalize[.03/(z + 2) Sqrt[z^2 + 2.5 z + 1.50 x]/((z - 1)^2 + .03)], x],
   z, Assumptions -> x > 0];

dzdx1 = ad /. {{z -> 0}, {z -> x}} // Differences // First // Simplify;

f2 = NDSolveValue[{
   f'[x] == 
     dzdx1 + (Rationalize[.03/(z + 2) Sqrt[z^2 + 2.5 z + 1.50 x]/((z - 1)^2 + .03)] /. z -> x), 
   f[0] == 0}, 
 f, {x, 0, 10}]

Replace f2[x][x] in functionsToPlot by f2[x] below in order to graph it.

Original

The first integral can be done once and for all with NDSolve. The second can be set up with ParametricNDSolve, which will reuse the functions it computes (one for each x). That will cut down on some time.

Clear[functionsToPlot, totalFunctions]
f1 = NDSolveValue[
    {f'[z] == 0.3/(z + 2) Sqrt[z^2 + 2.0 z + 0.7]/((z - 1)^2 + 1.1), 
    f[0] == 0}, f, {z, 0, 10}];
f2 = ParametricNDSolveValue[{
    f'[z] == .03/(z + 2) Sqrt[z^2 + 2.5 z + 1.50 x]/((z - 1)^2 + .03), 
    f[0] == 0}, f, {z, 0, 10}, {x}];
functionsToPlot[x_] := {
   (*1*) f1[x],
   (*2*) f2[x][x],
   (*3*) 0.5*Sqrt[x + 0.5]/x};

It's about ten times faster this way (0.58 sec. vs 5.46 sec.).

Note: One does have to pick an upper limit for x to use in NDSolve. I picked 10, but 3 is good enough for the example plot.

share|improve this answer
    
Nice. I was struggling with the second term when you posted this. –  Mr.Wizard Aug 19 at 21:02
    
@Mr.Wizard I was struggling with too, until just before I posted it. :) Thanks! –  Michael E2 Aug 19 at 21:04
    
This is very nice! That I have to pick an upper limit of x is no problem. I also have two cases (in the list) where the upper limit of the integral looks like (1-x)^2. Would that case also work with NDSolveValue or ParametricNDSolveValue? –  QuantumDot Aug 19 at 21:23
    
@QuantumDot If I understand you right, you can use NDSolveValue if the only x in the integral is in the limits -- it can be any function such as (1-x)^2. If x is in the integrand, too, then use ParametricNDSolveValue (unless the differentiate-integrate trick works). –  Michael E2 Aug 19 at 21:38
    
Yes, I just now figured out how it works (I had never heard of ParametricNDSolveValue before now). Thanks! –  QuantumDot Aug 19 at 21:43

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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