10
$\begingroup$

I'd like to be able to make InterpolatingFunctions for antiderivatives of functions that can't be integrated symbolically. However, the following code returns several error messages:

FunctionInterpolation[Integrate[Sqrt[1 + x^3], {x, 0, t}], {t, 0, 10}]

Here are the messages:

Thread::tdlen: "Objects of unequal length in {-1.25,-0.416667,0.416667,1.25}^{} cannot be combined. "

Thread::tdlen: "Objects of unequal length in {0.223144 +3.14159\ I,-0.875469+3.14159\ I,-0.875469,0.223144}\ {}\n cannot be combined. "

Thread::tdlen: "Objects of unequal length in {-1.25,-0.416667,0.416667,1.25}\ {} cannot be combined."

General::stop: "Further output of Thread::tlden will be suppressed during this calculation."

FunctionInterpolation::nreal: Near t = 1.25`, the function did not evaluate to a real number.

FunctionInterpolation::nreal: Near t = 1.3277777777777777`, the function did not evaluate to a real number.

What's going on? Is there a simple way to make this work? Changing Integrate to NIntegrate doesn't help, though the error messages are different:

NIntegrate::nlim: x = t is not a valid limit of integration.

NIntegrate::nlim: x = t is not a valid limit of integration.

NIntegrate::nlim: x = t is not a valid limit of integration.

General::stop: "Further output of NIntegrate::nlim will be suppressed during this calculation. "

$\endgroup$
  • 1
    $\begingroup$ your example is indeed integrated by MMA without problems an gives: (t^3)^(1/3)*Hypergeometric2F1[-(1/2), 1/3, 4/3, -t^3]. $\endgroup$ – Dr. Wolfgang Hintze Sep 17 '14 at 7:19
  • $\begingroup$ @Dr.WolfgangHintze Yes, it has no problems with the indefinite integral. The issue seems to be how I combined it with FunctionInterpolation. Strangely, other functions with more elementary integrals, e.g. Sqrt[1+x^2], do not produce any errors. $\endgroup$ – Jim Belk Sep 18 '14 at 1:28
  • $\begingroup$ (at)Jim: no, it's the definite integral. The indefinte integral is the true antiderivative, and it looks different. But, please have a look at my solution which shows (a) that in a (natural) 2-step procedure there are no problems with calculating the interpolation of your definite integral, and (b) this definite integral is different from the antiderivative. $\endgroup$ – Dr. Wolfgang Hintze Sep 19 '14 at 8:59
  • $\begingroup$ @JimBelk, This morning when I try first example, FunctionInterpolation[Integrate[Sqrt[1 + x^3], {x, 0, t}], {t, 0, 10}] does not give the errors you indicate, but FunctionInterpolation[Evaluate@Integrate[Sqrt[1 + x^3], {x, 0, t}], {t, 0, 10}] does give them. (The reason has to do with the ConditionalExpression returned when the Integrate is evaluated without assumptions or a definite value for t.) But can you confirm the same behavior? If not, what version of Mathematica are you using? $\endgroup$ – Michael E2 Sep 19 '14 at 10:35
  • $\begingroup$ (at) Michael E2: you might wish to look at my solution for Version 8: No problems encountered in interpolating both the definite integral (Jim's) and the antiderivative. $\endgroup$ – Dr. Wolfgang Hintze Sep 19 '14 at 19:33
11
$\begingroup$

Use NDSolve

antiD = NDSolveValue[{f'[x] == Sqrt[1 + x^3], f[0] == 0}, f, {x, 0, 10}]

Example usage:

Plot[antiD[x], {x, 0, 10}]

Mathematica graphics


Alternatively...

This works because this function can be antidifferentiated (by Mathematica).

antiD = FunctionInterpolation[
 Evaluate @ Integrate[Sqrt[1 + x^3], {x, 0, t}, Assumptions -> 0 < t < 10],
 {t, 0, 10}]

or...

integral[t_?NumericQ] := NIntegrate[Sqrt[1 + x^3], {x, 0, t}];
FunctionInterpolation[integral[t], {t, 0, 10}]

FunctionInterpolation evaluated its argument on a symbolic t. The pattern test ?NumericQ prevents evaluation of NIntegrate until an actual number is substituted for t. See also What are the most common pitfalls awaiting new users?. Note also that this way does many evaluations of NIntegrate, whereas the NDSolve method does just one integration.

$\endgroup$
  • $\begingroup$ That doesn't really answer the question. Why doesn't the NIntegrate form do the right thing (by substituting the appropriate numerical value for $t?$) $\endgroup$ – Igor Rivin Sep 17 '14 at 1:02
  • $\begingroup$ @IgorRivin I was adding that, but the OP didn't actually ask that Q. He asked about Integrate, and how to achieve his objective, which last Q I did answer. The first Q is probably not worth going into very deeply. Evaluate gets rid of the problem. $\endgroup$ – Michael E2 Sep 17 '14 at 1:19
  • 1
    $\begingroup$ When I need to generate accurate interpolation of a function I usually use Plot or Legacy'Plot for generating the points for Interpolation. FunctionInterpolation does not make adaptive refinement. $\endgroup$ – Alexey Popkov Sep 17 '14 at 4:21
  • 1
    $\begingroup$ @AlexeyPopkov That's a good point and I have done that, too. But NDSolve has PrecisionGoal and AccuracyGoal. I don't believe Plot has such controls (correct me, if I'm wrong -- I know it has "Refinement" and so forth, but it's a different way to control). The above and ifn = NDSolve[{y[t] == f[x[t]], x'[t] == 1, x[0] == 0}, y, {t, 0, 10}] are a pretty simple commands, and they tend to be the way I've done things the last few years. $\endgroup$ – Michael E2 Sep 17 '14 at 10:37
  • 1
    $\begingroup$ @AlexeyPopkov It adapts step size according to error estimates. The estimates can be fooled (a little), so it's not perfect. NDSolve sampling can be quite different from Plot. For instance NDSolveValue[{y[t] == Sin[x[t]], x'[t] == 1, x[0] == 0}, y, {t, 0, 2 Pi}] produces an interpolation with a max error of 6*10^-8 with 162 sample points whereas interpolating the Line produced by the corresponding Plot has a max error of 4*10^-6 with 431 sample points. But I don't want to exaggerate the differences; they're rather comparable, imo. $\endgroup$ – Michael E2 Sep 17 '14 at 13:53
4
$\begingroup$

There are two aspects in this problem depending on the exact definition of the task, and, as we shall see, both are completely solved by MMA without any additional facility.

The aspects are

a) calculate the interpolation of the definite integral $f=\int_0^t \sqrt{1+x^3} \, dx$

b) calculate the interpolation of the antiderivative $fad=\int\sqrt{1+x^3} \, dx$

Part a) interpolation of Jim's definite integral

There are no problems in interpolating your definite integral if you do it in two steps.

Step 1:

Calculate the definite integral (Mathematica does this exactly)

f = Integrate[Sqrt[1 + x^3], {x, 0, t}, Assumptions -> t >= -1]

(* -> (t^3)^(1/3) Hypergeometric2F1[-(1/2), 1/3, 4/3, -t^3] *)

Step 2:

Do the interpolation

fi = FunctionInterpolation[f, {t, 0, 10}]

(* -> InterpolatingFunction[{{0.`,10.`}},"<>"] *)   

And finally check the equality of f and fi:

Plot[{f, fi[t]}, {t, 0, 10}]

(* skip the picture *)

For completeness:

$Version

(* -> "8.0 for Microsoft Windows (64-bit) (October 7, 2011)" *)

This procedure works fine also with many other functions, e.g. integer powers other than 3 under the sqrt.

Part b) interpolation of the antiderivative

The antiderivative is given by the indefinite integral :

fad = Integrate[Sqrt[1 + x^3], x]

$\frac{2 \left(x+x^4+(-1)^{1/6} 3^{3/4} \sqrt{-(-1)^{1/6} \left((-1)^{2/3}+x\right)} \sqrt{1+(-1)^{1/3} x+(-1)^{2/3} x^2} \text{EllipticF}\left[\text{ArcSin}\left[\frac{\sqrt{-(-1)^{5/6} (1+x)}}{3^{1/4}}\right],(-1)^{1/3}\right]\right)}{5 \sqrt{1+x^3}}$

This quantity, besides having a small imaginary part

fad /. x -> 1.25

(* -> 2.34276 - 1.17815*10^-16 I *)

has a dicontinuity in both Re and Im:

Plot[{Re[fad], Im[fad]}, {x, 0, 3}, PlotRange -> {-2, 5}, 
 PlotLabel -> "Re and Im of the antiderivative fad of \!\(\*SqrtBox[\(1 + \
x^3\)]\)", AxesLabel -> {"x", "Re/Im[fad]"}]

enter image description here

Discontinuities of antiderivatives are quite common, and they have been discussed in this community several times (cf. for example Mismatch between numerical and analytic evaluation of an integral)

This could potentially cause a problem for FunctionInterpolate over an interval including the jump.

But it is not the case here !

Rather, the interpolation returns the correct complex valued function, with just announcing - as it should do - an accuracy problem connected with the jump:

fadi = FunctionInterpolation[fad, {x, 0, 5}]

FunctionInterpolation::ncvb: FunctionInterpolation failed to meet the prescribed accuracy and precision goals after 6 recursive bisections near x = {1.99219}. Continuing to refine elsewhere. >>

(* InterpolatingFunction[{{0.`,5.`}},"<>"] *)

Visualizing both Re and Im of the interpolation

Plot[{Re[fadi[t]], Im[fadi[t]]}, {t, 0, 3}, PlotRange -> {-2, 5}, PlotLabel -> 
  "Re and Im of the interpolation of the \nantiderivative fadi of \\!\(\*SqrtBox[\(1 + x^3\)]\)", AxesLabel -> {"x", "Re/Im[fadi]"}]

enter image description here

It shows agreement with the original function except, as expected, in the vicinity of the jumps.

Summarizing part b)

the interpolation of the true antiderivative is done perfectly well by Mathematica with only a justified accuracy error message due to the jump.

Conclusion

I have checked the complete code with a fresh kernel and can state that no problems are encountered in version 8 for both interpretations of the problem.

Best regards, Wolfgang

$\endgroup$

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.