# Tag Info

25

The simplest way to make new NIntegrate algorithms is by user defined integration rules. Below are given examples using a simple rule (the Simpson rule) and how NIntegrate's framework can utilize the new rule implementations with its algorithms. (Adaptive, symbolic processing, and singularity handling algorithms are seamlessly applied.) Basic 1D rule ...

19

Motivation (for a new semi-symbolic integration strategy) Consider the following integral, which cannot be done neigther by Integrate: Integrate[BesselJ[y, x^3], {x, 0, ∞}, {y, 0, 1}] (* Integrate[If[Re[y] > -(1/3), Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]), Integrate[BesselJ[y, x^3], {x, 0, Infinity}, Assumptions -> Re[y] <= -(1/3)]...

7

Because NIntegrate evaluates the integrands before starting the actual integration, in some cases (like this one) it is better to define the integrand function F with the signature F[x_?NumericQ]. BF[n_?NumericQ, x_?NumericQ] := BesselJ[n, x] NIntegrate[BF[9/2, x], {x, 0, 1}] (* 0.000148473 *) Integrate[BesselJ[9/2, x], {x, 0, 1}] %% // N (* Sqrt[2/\[...

6

Can NIntegrate remember or make full use of the result of a smaller upper-limit integral? Or generally, how to speed up the plot involving NIntegrate or is there any principle to do it? Below is a solution in the spirit of this request. In short, we find adaptively sampled points, compute integral estimates over the intervals having those points as ...

5

Actually the observed behavior is in full accord with the HoldAll attribute, just check what happens when there is no such attribute: nIntegrate[x + x, {x, 1, 2}] // Trace {{x + x, 2 x}, nIntegrate[2 x, {x, 1, 2}]} From the above it is seen that the arguments are evaluated before applying the rules associated with the function nIntegrate. The purpose ...

5

From the "Details and Options" section in the docs: NIntegrate first localizes the values of all variables, then evaluates f with the variables being symbolic, and then repeatedly evaluates the result numerically. So I guess this is expected behavior.

5

Here is a workaround suggested by the response I received from WRI: {sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., WhenEvent[Norm[{x'[t]}] < 0.0001, {x[t] -> 1.}]}, {x}, {t, 0, 200}, Method -> {"EquationSimplification" -> "Residual"} ] Plot[x[t] /. sol, {t, 0, 200}] Warning: This option works by converting the system to a DAE, ...

5

A bit speculative, but I think we can see why 171 is a magic number, if we factor the Gamma[n] from the integral: Table[ {n, NIntegrate[Gamma[x] /Gamma[x + n], {x, 1, 2}]} , {n, 168, 173}] {{168, 7.20924*10^-304}, {169, 4.26119*10^-306}, {170, 2.41873*10^-308}, {171, 0.}, {172, 0.}, {173, 0.}} 171 is where the integral fails as you can see the ...

5

What is observed for testfunc and integrand2 is explained with the use of adaptive sampling, symbolic processing, and singularity handlers by NIntegrate. The home-cooked Simpson integration strategy in the question is too simple for these integrands. For the function integrand it seems that only the adaptive sampling gives the advantage of NIntegrate. Needs[...

5

mat[x_?NumberQ] := {{x^2, 1, 0}, {1, x^2, -1}, {0, -1, x}} ei[x_] := Eigensystem[mat[x], 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}][[1, 1]] NIntegrate[ei[x], {x, 0, 4}] (* 8. *)

5

First we shall define 'the integration on a curve'. Traditionally, this is defined as integration of f.dl where dl is the length of a small part of the curve. So, using t as a medium, we can explicitly write out the curve's function on a complex plane, here let's assume it's z=2 Exp[I t]. Then we can use t, a real number, as the integration variable, which ...

5

Maybe commandeer FindRoot with the Villegas-Gayley trick: Updated, with the order of the steps taken by FindRoot saved in icsteps. The results of FindRoot, as saved by NDSolve and shown below as DownValues[], have been sorted by Mathematica and are not in the order in which there were called. This update stores the order in icsteps. Clear[x, t]; Internal...

4

There is some confusion using Plot. Compare the plot of the non-evaluated y in the question with the plot of the evaluated y: Grid[{{Plot[y[0.25, w - 2], {w, 0, 4}, PlotRange -> All, ImageSize -> Medium], Plot[Evaluate@y[0.25, w - 2], {w, 0, 4}, PlotRange -> All, ImageSize -> Medium]}}] We can see that going through the ...

4

Here is another type of integration rule. It has a customized error estimator, incorporating which was my main interest. In fact, the rule below computes an integral of arbitrary dimension, using its own internal formula for the integral as well. Introduction A rule of the type NIintegrateGeneralRule[{abscissas, weights, errweights}] computes the ...

4

I think this gives what you want: We construct the integral inside the definition of H by adding another ODE to the NDSolve system, which I called logH. This in fact calculates the integral from ic, not from 0. So to define H we need to subtract logH[0] from logH[t] before exponentiating. This should be a much more accurate (and faster) way of computing ...

3

You have many superfluous sets of {} that generate unexpected output in your code. In particular, the Interpolation functions generated by NDSolve were not returning a scalar value, but instead a unidimensional vector, i.e. a list containing a single value instead. That was probably an unintended consequence of the extra sets of braces in the definitions of ...

3

Clear["Global`*"] yy = 10^-4; rr = 0.999; xx = 10^-15; zz = 10^-4; mm = 10^-4; yy + rr + xx^2 + zz - mm^2 - zz^2/24 ic = -17.5 s = NDSolve[{D[y[t], t] == (3 y[t])/5 - (12 m[t]^2 y[t])/5 + (2 r[t] y[t])/ 5 - (6 x[t]^2 y[t])/5 + (3 y[t]^2)/5 + (7 y[t] z[t])/ 5 - (y[t] z[t]^2)/10, D[r[t], t] == -((2 r[t])/5) - (12 m[t]^2 r[t])/5 + (2 r[...

3

As I mentioned in comments, your main problem was the placement of the definition of z; that definition does not need to be re-evaluated every time there is a modification within the body of Manipulate since it's already delayed. The correct placement for such a definition would be in the Initialization code for the Manipulate. Similarly, the definitions of ...

3

Try: {sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., WhenEvent[Norm[{x'[t]}] < 0.0001, x[t] -> 1.; x[t] -> 1]}, {x}, {t, 0, 200}]

3

Fixing it ... Is there any settings that I can turn on to fix this? Does using larger working precision produce results you expect? F[n_] := NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}, WorkingPrecision -> 60, MinRecursion -> 6, MaxRecursion -> 20] In[10]:= F[170] Out[10]= 0....

3

This answer is just a response to the following part of the question: ……generally, how to speed up the plot involving NIntegrate or is there any principle to do it? When the symbolic processing of NIntegrate isn't necessary for getting the correct result of integration, "SymbolicProcessing" -> 0 is also a general way for speeding up: f[t_] := ...

3

Expressing the entire business in terms of SphericalBesselJ[] cures the problem: NIntegrate[r With[{x = r BesselJZero[15/2, 1]}, Sqrt[2 x/π] SphericalBesselJ[7, x]]^2, {r, 0, 1/50}] 1.1879560281974252*^-27 The nice thing about SphericalBesselJ[] is that it does not auto-evaluate to a potentially numerically unstable combination of ...

3

int[a_, b_] = Integrate[Sin[y1[t] y2[t] r]/r, {r, a, b}, Assumptions -> a < b] odes = {y1'[t] == y2[t] y3[t], y2'[t] == -y1[t] y3[t] + int[y2[t], y3[t]], y3'[t] == -0.51 y1[t] y2[t] t, y1[0] == 0, y2[0] == 1, y3[0] == 2} sol = NDSolve[odes, {y1, y2, y3}, {t, 0, 12}] Edit to answer your question in your comment: int[a_?NumericQ, b_?NumericQ,...

2

I tried to make it a comment and then can't control the words. I think the problem is arising from the fact that the value of the function and its derivative is too small to near x=0. Plot[BesselJ[9/2, x], {x, 0, .01}] Plot[Evaluate[D[BesselJ[9/2, x], x]], {x, 0, 0.01}] As you can see the derivative is highly oscillatory at this small range. When ...

2

EDIT It can get more complicated. Look at a contour given by r[t_] := 1 + 2 Cos[t]; Original post The result of the integral for a (reasonable) closed contour is just 2 I times the area enclosed by the contour. Proof: the integral is $$\int \left(z^*+z\right) \, dz$$ Letting $$z=x+i y$$, $$dz=dx+i dy$$ the integral becomes \int 2 x (dx + i dy)=...

2

The DE has problems at y[x] == 0 as noted, but also at x == 1, where there is a pole of order 2. This suggests that numerical issues near the singularities could cause substantial error to accumulate. One can use Piecewise to substitute the limiting value at a discontinuity, but this won't cure any numerical instability in the neighborhood of the ...

2

As requested, comment made into an answer. The problem was with the integration limits in the OP. The points on the circle satisfy (x-x0)^2 + (y-y0)^2 == R^2 so for a point x in the range [x0-R,x0+R] the limits on y are {y0 - Surd[R^2 - (x-x0)^2,2],y0 + Surd[R^2-(x-x0)^2,2]} I use Surd here because MMA may simplify more easily than with Sqrt, since ...

2

Let us have a single integrand function first: F[x_?NumericQ, xp_?NumericQ] := ( EllipticK[(4 x xp)/(x + xp)^2] Sinh[x] Sinh[ xp])/((x + xp) (Cosh[202] + Cosh[2 x]) (Cosh[202] + Cosh[2 xp])); obtained by expanding f[x, 101.] f[xp, 101.] g[x, xp]. For me (with Mathematica v. 10.3), the integral fails with NIntegrate::inumri: "The integrand (....

2

This works Clear[x, y, z] sol = NDSolve[Join[ Thread[{x'[t], y'[t], z'[t]} == {-z[t] x[t] + 1.5 Log[z[t]], z[t] - 1, -(0.64) y[t] + (x[t] + Log[z[t]])/0.8}], Thread[{x[0], y[0], z[0]} == {1, 1, 1}]], {x, y, z}, {t, 0, 100}]; ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 100}]

2

I have revised the code to use StreamPlot instead of VectorPlot and NDSolve. The latter command created issues when the domain for the solution extended to negative times. I don't know why, but NDSolve produced integral curves that did not follow the direction field produced by VectorPLot as provided by MarcoB's nice reformulation of my original code. I will ...

Only top voted, non community-wiki answers of a minimum length are eligible