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!

Improve this question
$\endgroup$
4
  • 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$ Commented May 8, 2019 at 17:34
  • 1
    $\begingroup$ Search the Mathematica documentation for 'NIntegrate'. It is easy to find. $\endgroup$ Commented May 8, 2019 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$ Commented May 8, 2019 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$ Commented May 10, 2019 at 3:06

3 Answers 3

Reset to default
5
$\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

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ You could also interpolate your table to produce a function. $\endgroup$ Commented May 8, 2019 at 18:19
  • 1
    $\begingroup$ @Michael E2 You're right. I should have done it right away! $\endgroup$ Commented May 8, 2019 at 20:49
  • 1
    $\begingroup$ Recommend that you use the option DataRange -> {0, 1} with the ListLinePlot. $\endgroup$ Commented May 10, 2019 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$ Commented May 12, 2019 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

Improve this answer
$\endgroup$
2
$\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

Improve this answer
$\endgroup$
2
  • $\begingroup$ Related to the matter of integrating a Chebyshev series, have you seen this already? $\endgroup$ Commented May 12, 2020 at 11:57
  • $\begingroup$ @J.M. No, it seems I missed that. (+1) I'm familiar with the idea, and I've read the El-Gendi paper the OP cites. $\endgroup$ Commented May 12, 2020 at 13:39

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.