2
$\begingroup$

The function (in raw input form) and bounds I used are as follows:

Integrate[(0.3950832348257582*Sqrt[(-(-1 + z))*z]*
(-1.8816764231589205 - 15.31803072355397*z + 55.36247428645651*z^2 - 
57.24415070961543*z^3 + 19.08138356987181*z^4 + 
(13.642154067902172 - 8.202924565932532*z - 43.60199664171326*z^2 + 
57.24415070961543*z^3 - 19.08138356987181*z^4)/
E^(1.7146776406035664/(1.*z - 1.*z^2))))/
E^(0.27434842249657065/(1.*z - 1.*z^2))/((1. - 1.*z)^2*(-1. + z)*z), 
{z, 0, u}]

The function was originally a ConditionalExpression, which I threw into the Normal[] function in hopes of making it solvable. The output that was returned upon trying to evaluate the integral is identical to the input (the above snippet).

u is an arbitrary value between 0 and 1 that will be eliminated in another integral later in my calculations.

Any help is much appreciated!

$\endgroup$
  • 2
    $\begingroup$ People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this meta Q&A helpful $\endgroup$ – Michael E2 May 8 at 17:34
  • 1
    $\begingroup$ Search the Mathematica documentation for 'NIntegrate'. It is easy to find. $\endgroup$ – rmw May 8 at 17:37
  • 2
    $\begingroup$ Integrate is an exact solver, and sometimes round-off error from floating point coefficients make it fail. In this case, even if you Rationalize[] the coefficients, Integrate fails. $\endgroup$ – Michael E2 May 8 at 18:22
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 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! $\endgroup$ – Michael E2 May 10 at 3:06
4
$\begingroup$

Maybe this helps?

nsol = Table[
  NIntegrate[(0.3950832348257582*
       Sqrt[(-(-1 + z))*z]*(-1.8816764231589205 - 
         15.31803072355397*z + 55.36247428645651*z^2 - 
         57.24415070961543*z^3 + 
         19.08138356987181*
          z^4 + (13.642154067902172 - 8.202924565932532*z - 
            43.60199664171326*z^2 + 57.24415070961543*z^3 - 
            19.08138356987181*z^4)/
          E^(1.7146776406035664/(1.*z - 1.*z^2))))/
     E^(0.27434842249657065/(1.*z - 1.*z^2))/((1. - 1.*z)^2*(-1. + z)*
      z), {z, 0, u}], {u, 0, 1, 0.01}]

ListLinePlot[nsol, GridLines -> Automatic]

enter image description here

f[u_] = Interpolation[Thread@{Table[u, {u, 0, 1, 0.01}], nsol}, u];
Plot[f[u], {u, 0, 1}, GridLines -> Automatic]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ You could also interpolate your table to produce a function. $\endgroup$ – Michael E2 May 8 at 18:19
  • 1
    $\begingroup$ @Michael E2 You're right. I should have done it right away! $\endgroup$ – rmw May 8 at 20:49
  • 1
    $\begingroup$ Recommend that you use the option DataRange -> {0, 1} with the ListLinePlot. $\endgroup$ – Bob Hanlon May 10 at 22:13
  • $\begingroup$ @Bob Hanlon Thanks for the hint. I could have chosen f = ListInterpolation[nsol, {{0, 1}}] would have been faster. Thanks again! $\endgroup$ – rmw May 12 at 9:53
4
$\begingroup$

You could use NDSolveValue:

yp[z_?NumericQ] := 
  Piecewise[{{(0.3950832348257582*
         Sqrt[(-(-1 + z))*z]*(-1.8816764231589205 - 
           15.31803072355397*z + 55.36247428645651*z^2 - 
           57.24415070961543*z^3 + 
           19.08138356987181*
            z^4 + (13.642154067902172 - 8.202924565932532*z - 
              43.60199664171326*z^2 + 57.24415070961543*z^3 - 
              19.08138356987181*z^4)/
            E^(1.7146776406035664/(1.*z - 1.*z^2))))/
       E^(0.27434842249657065/(1.*z - 1.*z^2))/((1. - 1.*z)^2*(-1. + 
          z)*z), 0 < z < 1}}
   ];
integral = NDSolveValue[{y'[z] == yp[z], y[0] == 0}, y, {z, 0, 1}]

Plot[integral[u], {u, 0, 1}]

enter image description here

$\endgroup$
1
$\begingroup$

If you want a highly accurate approximant, a Chebyshev approximation is a good approach. See Trefethen, Approximation Theory and Approximation Practice, Boyd, Solving Transcendental Equations, and this answer by J.M; a Chebyshev series may be antidifferentiated with iCheb.

Here is the basic approximation. Since Chebyshev polynomials stay between $\pm1$, the error can be estimated from the coefficients of the tail of a rapidly convergent Chebyshev series. The plot shows when the coefficients run into the round-off error limit, which is around machine epsilon times the maximum absolute coefficient. The horizontal gridline at the bottom of the plot shows the smallest error one could hope for; round-off error tends to be somewhat larger than this.

yp[z_?NumericQ] := (* OP's function with discontinuities at 0,1 removed *)
  Piecewise[{{(0.3950832348257582*
         Sqrt[(-(-1 + z))*z]*(-1.8816764231589205 - 
           15.31803072355397*z + 55.36247428645651*z^2 - 
           57.24415070961543*z^3 + 
           19.08138356987181*
            z^4 + (13.642154067902172 - 8.202924565932532*z - 
              43.60199664171326*z^2 + 57.24415070961543*z^3 - 
              19.08138356987181*z^4)/
            E^(1.7146776406035664/(1.*z - 1.*z^2))))/
       E^(0.27434842249657065/(1.*z - 1.*z^2))/((1. - 1.*z)^2*(-1. + 
          z)*z), 0 < z < 1}},
   0
   ];

deg = 256;
chebnodes = N[Rescale[Sin[Pi/2 Range[-deg, deg, 2]/deg]]];
yvals = yp /@ chebnodes // Quiet;
chebcoeffs = Sqrt[2/deg] FourierDCT[yvals, 1];
chebcoeffs[[{1, -1}]] /= 2;
ListPlot[RealExponent[chebcoeffs], 
 GridLines -> {None, {Max@Abs@chebcoeffs*$MachineEpsilon // RealExponent}}]

Below is the iCheb routine from the linked answer above, which computes the Chebyshev series of an antiderivative of a given series. The constant of integration needs to be computed from the initial antiderivative, and the first Chebyshev coefficient needs to be adjusted accordingly. We can trim the coefficients of the tail that are below round-off error. This step is optional and make computing with the Cheybshev series only slightly more efficient.

(*Integrate a Chebyshev series-- cf.Clenshaw-Norton,Comp.J.,1963,p89,eq.(12)*)
Clear[iCheb];
iCheb::usage = "iCheb[c, {a, b}, k] integrates the Chebyshev series c, plus k";
iCheb[c0_, {a_, b_}, k_: 0] := Module[{c, i, i0}, c[1] = 2 First[c0];
  c[n_] /; 1 < n <= Length[c0] := c0[[n]];
  c[_] := 0;
  i = 1/2 (b - a) Table[(c[n - 1] - c[n + 1])/(2 (n - 1)), {n, 2, Length[c0] + 1}];
  i0 = i[[2 ;; All ;; 2]];
  Prepend[i, k - Sum[(-1)^n*i0[[n]], {n, Length[i0]}]]]

ClearAll[trimCC];
trimCC[cc_] := 
  With[{drop = 
     1 - With[{m = Max@Abs@cc}, 
       Module[{err = 0.}, 
        LengthWhile[
         Reverse@cc, (err += Abs[#]) < $MachineEpsilon*m &]]]},
   Drop[cc, -drop] /; drop > 2];
trimCC[cc_] := cc;

intcc = iCheb[chebcoeffs, {0, 1}];
intcc[[1]] += intcc.(-1)^Range[Length@intcc];  (* adjust constant of integration *)
intcc = trimCC[intcc];
intCS[u_] := intcc.Cos[Range[0, Length@intcc - 1] ArcCos[2 u - 1]];

Plot[intCS[u], {u, 0, 1}]

enter image description here

In comparison with directing numerical integration, which is relatively slow (we raise the PrecisionGoal slightly to get a more accurate numerical integral), it's a pretty good approximation!:

Plot[intCS[u] - 
  NIntegrate[(0.3950832348257582*
       Sqrt[(-(-1 + z))*z]*(-1.8816764231589205 - 
         15.31803072355397*z + 55.36247428645651*z^2 - 
         57.24415070961543*z^3 + 
         19.08138356987181*
          z^4 + (13.642154067902172 - 8.202924565932532*z - 
            43.60199664171326*z^2 + 57.24415070961543*z^3 - 
            19.08138356987181*z^4)/
          E^(1.7146776406035664/(1.*z - 1.*z^2))))/
     E^(0.27434842249657065/(1.*z - 1.*z^2))/((1. - 1.*z)^2*(-1. + z)*
      z),
   {z, 0, u}, PrecisionGoal -> 12, AccuracyGoal -> 16], {u, 0, 1}]

enter image description here

$\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.