# Tag Info

## New answers tagged numerical-integration

1

I think substituting a in the definition of int, which is called by a is trouble, and what you really need is to make sure int isn't evaluated when z[t] isn't a number. The following seems to work, albeit not very quickly: Clear[a]; f[t_] := Sin[t^2]; int[t_?NumericQ, z_?NumericQ] := NIntegrate[Exp[-f[t]*z*v^2], {v, 0, 1}]; a = NDSolve[{ x'[t] == -3*(x[...

0

Somewhat similar to J.M.'s but using Bezier curves instead of B-splines. simpsonSegment[{{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}] := FilledCurve[ {Line[{{x1, 0}, {x1, y1}}], BezierCurve[{{x1, y1}, {x2, 1/2 (-y1 + 4 y2 - y3)}, {x3, y3}}], Line[{{x3, y3}, {x3, 0}}]} ]; Manipulate[ Plot[f, {x, a, b}, PlotStyle -> {Red, Thick}, AxesOrigin -&...

3

This evaluates without errors: z[x_] := 2458.31 - 100.087 x + 1.23213 x^2 - 0.0046743 x^3 lamavg[t_] := Min[1, 0.01 + 0.07 z[t]] Ufit2[M_, t_] := 0.50519 + 3.127*10^10/M^2 - 274337/M + 2.12127*10^-10 M - 1.92858*10^-20 M^2 - 6.20762*10^-11 t; e = 1/100; l = (126/100)*10^(31); sol = 3*^8; DifEq = D[P[M, t],t] == -M l/ sol^2 D[(1 - e)/e ...

0

I encountered a problem copying and pasting your code, as HypergeometricU and LaguerreL begin on new lines, but I believe you meant for them to be part of the definitions for DI1 and DI2. Fixing that, the problem of non-numerical integrand comes from undefined global U, which can be set inside Block, and FindRoot attempting to evaluate its argument ...

1

You can find and exclude the part of the domain where the integrand is discontinuous. To find where the integrand is discontinuous, look for where the argument of Arg is real and negative, since Arg has a branch cut discontinuity there: z[\[Gamma]_,\[Beta]_,\[Omega]_]:=(Times@@(secondOrderDelayedAllPass[#,\[Omega]]&/@\[Gamma]))* ...

2

This can be written down as the parametric iterated integral $$I(a)=\int_0^a\left( \int_0^x f(y)\mathrm dy\right) \mathrm dx$$ NIntegrate supports the syntax NIntegrate[f,{x,xm,xM},{y,ym,yM}] (where $x_m,x_M$ are the bounds of the outer integral, $y_m,y_M$ of the inner one) to compute this kind of integrals. This means that you can evaluate your result ...

2

This is not a complete answer, the code provides results that demonstrate the need for more detailed investigation. Re-definition The re-defintion uses exact numbers and adds options argument to F. Clear["Global*"]; χ[x_] := (1 - x^2) Exp[-x^2/2] v1[x_, T_, A_, ν_] := Sqrt[2 ν A] (-Tanh[Sqrt[A/2 ν] (x - Sqrt[-Sqrt[A] T])] + Tanh[Sqrt[A/2 ν] (x + ...

3

Well, since an example NDEigensystem was not provided, I went ahead an took one from the docs. {vals, funs} = NDEigensystem[ {-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u, {x, y} ∈ RegionDifference[Cuboid[{-3, -3/2}, {0, 3/2}], Disk[]], 10, Method -> {"SpatialDiscretization" -> {"FiniteElement", {"...

2

I think you want something like this. I am assuming in your actual application f will require numerical integration, even though in this example it doesn't. f[x_, y_] := x^2 y + x y^2 g[y_?NumericQ] := NIntegrate[f[x, y], {x, 0, 1}] Plot[g[y], {y, 0, 1}, AxesLabel -> {y, g[x]}]

2

The short answer is no, the result is insignificant (as is). The error, according to AccuracyGoal -> 10 is around 10^-10, which is 100 times larger than the result. To check 10^-12, you should bump AccuracyGoal up above that, say to 12 plus half MachinePrecision, or 12 + 8. NIntegrate[ E^(-ω/(2 a)) E^(-w/(2 b)) ((1 - Cos[(ω + w + Ω) t]) / (Sqrt[2] (...

4

This answer shows how to define a new NIntegrate rule that evaluates f in the list of two integrands {f[x],g[f[x]]+h[x]} only once per sampling point. The answer can be also easily modified into an answer of "NIntegrate over a list of functions". The definition of the NIntegrate rule LessEvaluationsRule given below is also aimed to be didactic and ...

4

First, NIntegrate[f1[x], {x, xmin, xmax}] usually proceeds by constructing an ExperimentalNumericalFunction from the expression for f1[x]. This will circumvent an attempt to memoize f1 in the OP's manner, f1[x_] := f1[x] =.... One can prevent this by memoizing the function with ?NumericQ checks via f2[x_?NumericQ] := f2[x] = .... One thing to consider is ...

-1

Testing my comment does indicate that one can cut the number of evaluations of f[] in half easily. f[x_, y_, z_] := Module[{}, totCalls++; Exp[Sin[x]] + Cos[y + z] ] totCalls = 0; NIntegrate[ {f[x, y, z], Sqrt[f[x, y, z]] + x}, {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, ...

3

I got around to actually evaluating your code and I realized that g is not remembering its values; DownValues[g] only has a length of three. The "solution" is to restrict the function to numeric values, per The difference between "SymbolicProcessing" -> 0 and restricting the function definition to numeric values only, but doing that actually ...

1

Of course in this case you can trust the symbolic result of Integrate, but the point raised by the question becomes especially important when there is no analytical solution. Here is an example suffering from an even more rapid decay that causes the same numerical problems. Finding the right options for NIntegrate isn't so obvious: a = 10^6; integrand = ...

2

Most of the value of the integral occurs before $w\le 0.00002$. NIntegrate[Exp[-w/a] t, {w, 0, 0.00002}, WorkingPrecision -> 16] 0.009999999979388476 At MachinePrecision, the total value of the integral can be found at $w\le 0.1$. NIntegrate[Exp[-w/a] t , {w, 0, 0.1}] 0.01 It is therefore not just a question of MaxRecursion: NIntegrate[...

0

RegionPlot does indeed produce an error. I have contacted Wolfram support with this example. They suggested a workaround using Plot3D and RegionFunction. r1 = 0.00029938135; r3 = 2.507341314*^-8; r4 = 0.00001129433024; r5 = 0.005435396430; r8 = 0.5761520216; cs[t_?NumericQ, τ_?NumericQ, d_?NumericQ] := Exp[t - t*r5]* NIntegrate[Exp[ξ*(-1 + r5)]*(d*Exp[τ*...

2

I did manage to get a 6-7 times speed up by separating the definition of psi into two integrals and experimenting with the option values. (I am not sure how significant the speed-up is... OP mentioned in the comments that 1 second is the goal.) Here are the two new functions: Clear[psi1] psi1[x_?NumericQ, y_?NumericQ, z_?NumericQ, t_?NumericQ, opts : ...

4

There are several ways to handle singularities with NIntegrate. Concerning your question, Since you already know the location of singularities, simply remove them with Exclusions during integration. NIntegrate[E^(-ω/(2 a)) E^(-w/(2b))((1 - Cos[(ω +w + Ω) t ])/(Sqrt[2] (ω + w + Ω)(ω + Ω))) Sqrt[w ω],{w, 0,1},{ω, 0, 1}, MaxRecursion -> 300, AccuracyGoal -&...

4

TL;DR Use HeavisideTheta's properties before integration. This is my strategy. First the HeavisideTheta gives you the following integration limits: $$0\leq y \leq 1-x \qquad \& \qquad 0\leq x \leq 1$$ $$0\leq x \leq 1-y \qquad \& \qquad 0\leq y \leq 1$$ In both cases I used Integrate first then NIntegrate. In the first case I could not ...

5

The shortest and best way between two truths of the real domain often passes through the imaginary one. — Jacques Hadamard By taking a complex path, I get the answer without any complaints from Mathematica. parabolic[a_, x_] = Simplify[InterpolatingPolynomial[{{-10, 0}, {0, a}, {10, 0}}, x]] With[{a = 1}, Re[NIntegrate[With[{x = x + I ...

3

You can try setting the AccuracyGoal lower than WorkingPrecision which yields the correct result without a reported warning. NIntegrate[Exp[-x^2] Cos[100 x], {x, -10, 10}, Method -> {"LevinRule"}, WorkingPrecision -> 50, AccuracyGoal -> 35] 5.1113608752199120138254477520179596033660767259737*10^-46 NIntegrate[Exp[-x^2] Cos[100 x], {x, -10, ...

1

I cannot evaluate your code due to missing functions but here is a shot in the dark: CMB[a_?NumericQ, lambda_?NumericQ] := FindRoot[ns[i, lambda, a] - .96 == 0, {i, 15}] END[a_?NumericQ, lambda_?NumericQ] := FindRoot[eps[i, lambda, a] - 1 == 0, {i, 1}] Ne[lambda_?NumericQ, a_?NumericQ] := With[{min = END[a, lambda], max = CMB[a, lambda]}, ...

1

Table[NIntegrate[ im[s]/(s (s - SetPrecision[si[[i]], MachinePrecision])), {s, 4 m^2, 2}, Method -> "PrincipalValue", Exclusions -> si[[i]] == s, MaxRecursion -> 500], {i, 1, 21}] {2.70869, 0.196226, -0.196659, -0.300094, -0.328605, -0.331732, -0.325217, -0.315038, -0.303723, -0.292415, -0.281655, -0.271706, -0.262715, -0.254788, ...

5

You need to construct an event for each i from 1 to n: Block[{n = 2, a = 1.1}, vars = Table[x[i], {i, n}]; eqns = Table[x[i]'[t] == a - x[i][t], {i, n}]; initcond = Table[x[i][0] == 0.3*i, {i, n}]; evts = Table[With[{i = i}, WhenEvent[x[i][t] == 1, x[i][t] -> 0]], {i, n}]; sol = NDSolve[{eqns, initcond, evts}, vars, {t, 0, 10}]; ] Plot[Evaluate[...

5

Aside from getting around the apparent weakness in the "LevinRule"* as others have suggested, here is another way to verify the total probability is 1, namely, by changing variables. {transformation} = Solve[{u1, u2} == {Log[x1/(1 - x1 - x2)], Log[x2/(1 - x1 - x2)]}, {x1, x2}, Reals] (* {{x1 -> E^u1/(1 + E^u1 + E^u2), x2 -> -E^-u1 (-E^u1 + E^...

1

Just picking up the idea mentioned in the comments, your code would be f[r_] := 1/(1 - 4/(r*Sqrt[Pi])*Gamma[1.5, r^2/(4*0.2)]) F[y_?NumericQ] := NIntegrate[f[r], {r, 1, y}] LogLogPlot[F[r], {r, 1, 100}]

2

Here is code that makes the plot of 1st ReandIm of the function without messages. Clear[f, "G*", ϕ] f[x_?NumberQ] := Exp[-(x - 5)^2] G1[b_?NumberQ, σ_, λ_] := 0 G2[b_?NumberQ, σ_, λ_] := (1/π) Sqrt[ b/σ] EllipticK[ Abs[(λ^2 - 4 (σ - b)^2)/(16 σ*b)]]; G3[b_?NumberQ, σ_, λ_] := (4/π)*((b)/(Sqrt[\ λ^2 - 4 (σ - b)^2])) EllipticK[ Abs[(16 σ*b)/(λ^2 ...

1

Not a real answer, too long to put in a comment... What I would suggest is to use smaller precision goals and restrain NIntegrate's adaptive algorithm. This way you can get an answe in a relatively short time, analyze does it make sense, and continue with refinements that require more computational time. The code below produces results under 35s on my ...

12

The sampling points are insufficient in the first rule applications. Increasing them, say, with MinRecursion produces the expected result: NIntegrate[f[x1, x2, 0, 0, 1, 1, 0], {x1, 0, 1}, {x2, 0, 1 - x1}, MinRecursion -> 1] (* 1. *) An alternative is to use Cartesian rules: NIntegrate[ f[x1, x2, 0, 0, 1, 1, 0], {x1, 0, 1}, {x2, 0, 1 - x1}, Method ...

2

Having a helper function rhs, which evaluates only with a numeric vector as argument, for the right-hand side of the force equation lets you use vectors as you want. This way the undesired symbolic precalculation (threading of drag (v.v) Normalize[v] with {0, 0, gravity}) is bypassed and the solving continues numerically. See this answer for a bit more ...

2

The interior integral, done by NDSolve. We save the results for speed. Clear[i1]; i1[u_?NumericQ] := i1[u] = NDSolveValue[{ff'[s] == 2 q[s, u]/d[s, u], ff[0] == 0}, ff, {s, -0.28, 0.28}] The outside integral "SymbolicProcessing" -> 0 is for speed, which it seems to give in this case, and the interval for u is assumed to be ±0.1, as it was for x in ...

2

ClearAll["Global*"] Remove["Global*"] G = 1/100; ωc = 15; f[t_ , ω_] := G ω Exp[-ω/ωc] (Sin[(ω - 1) t/2]/((ω - 1) t/2))^2 Int[t_?NumericQ] := NIntegrate[f[t , ω], {ω, 0, Infinity}] Plot[1 - 2*t^2*Int[t], {t, 10^-15, 100}] MMA has problems with function Sinc give some error messages: f[t_ , ω_] := G ω Exp[-ω/ωc] (Sinc[(ω - 1) t/2])^2 Int2[t_?...

1

Here is a way to do it: (* Needs["NDSolveFEM`"]; \[CapitalOmega]=Rectangle[{-0.05,-0.1},{0.2,0.1}]; mesh=ToElementMesh[\[CapitalOmega],"IncludePoints"\[Rule]{{0,0}}]; *) sol = NDSolveValue[{-D[u[x, y], x, x] - D[u[x, y], y, y] - u[x, y] == 0, DirichletCondition[u[x, y] == 1, x == 0 && y == 0], DirichletCondition[u[x, y] == 1, y == ...

0

Sinc^2[_] is wrong and should be Sinc[_]^2 unless you mean Sinc[Sinc[_]]. integral1 should be defined as a function because it depends on t and ω. For some reason NIntegrate happen to transform the integrand through ω into a non-numeric expression (which is obviously wrong) unless it is prevented by the syntax ω_?NumericQ. G = 0.01; ωc = 15; β = 1; ω0 = 11;...

1

One can use NDSolve to construct a particular "indefinite" integral, e.g., $$f(z) = \int_0^z f'(ζ) \; dζ$$ First define integrand to be the OP's expression. The the following computes its integral. Clear[df]; df[z0_?NumericQ] := df[z0] = NIntegrate[integrand /. z -> z0, {ϕ, 0, 7 Pi/18}, Method -> {"GaussKronrodRule", "Points" -> 11}, ...

1

Spaces matter. If we check the OP's integral4 we find in the output an undefined symbol βω. Probably it was meant to be a β ω, meaning β * ω. integral4 := G ω Exp[-ω/ωc] ((1 - Cos[ω τ])/ω^2) Coth[βω/2] integral4 (* (0.01 E^(-ω/15) (1 - Cos[τ ω]) Coth[βω/2])/ω *) integral4 = G ω Exp[-ω/ωc] ((1 - Cos[ω τ])/ω^(2)) Coth[β ω/2]; Plot[-(1/τ) Log[ 1/2 +...

10

Update-1: The initial conditions in the question were wrong/incomplete. (removed 1/0 errors) Update-2: The 1D Euler equations were modified to match this source. (ultimately not necessary in V10.4, but is in V8) Update-3: Method options in NDSolve were modified to produce an accurate result. (ENO schemes are not yet supported, but the proposed answer below ...

4

[...] I am only interested in very fast numerical methods, no analytical results are needed. [...] I have no idea how I can do it in Mathematica The package AdaptiveNumericalLebesgueIntegration.m has Lebesgue integration strategy and rules implementations and it is discussed in detail in the blog post "Adaptive numerical Lebesgue integration by set ...

0

The documentation of FindRoot says FindRoot first localizes the values of all variables, then evaluates f with the variables being symbolic, and then repeatedly evaluates the result numerically. So the system first evaluates Integrate[(x^{0, 2, 4})*f[x, a, b, c], {x, -Infinity, Infinity}] - {1, 2, 10} // N which gives the errors you see (...

2

The integrand cannot be solved when $|z|\rightarrow1$ Therefore, let's integrate the function on a possible domain, numerically. zdat=Table[NIntegrate[f, {Phi, 0, 7 Pi/18}, {z, 0, i}], {i, -0.95, 0.95, 0.1}]; This gives us a list of values, from which we can approximate a function for this part of the domain for z. Plotting this list gives us: lp = ...

16

Yes you can. Below is a fairly general, Mathematica-compiled, fast and robust version. Examples 1. Michaelis-Menten kinetics Michaelis-Menten kinetics for enzyme-directed substrate conversion. The enzyme (e) converts the susbtrate (s) through an enzyme-substrate complex (c) to the product (p). For comparison, I've included the deterministic ODE system ...

3

With your definitions, use: integral[v_?NumericQ, x_?NumericQ] := (1/d[v, x])*Exp[NIntegrate[2 q[s, x]/d[s, x], {s, 0, v}]] Plot3D[ integral[v, x], {v, -0.5, 0.5}, {x, -0.1, 0.1}, PlotPoints -> 10, MaxRecursion -> 0 ] Read this FAQ to see why you need to use NumericQ in this case. Having said that, your function assumes insanely high ...

1

There's not enough information in the question, but this approach fits one interpretation of the description of the situation. One can do the integral in the original NDSolve calls that generated rad and mm. (This may be what J.M. was suggesting in a comment. Or he may have meant executing a separate NDSolve call. I would recommend the approach given ...

0

If you solve for the second derivatives, you won't have to use "EquationSimplification" -> "Residual" and things will work ok. Solving for the second derivatives be faster if you start with exact coefficients. Also, if you clear l, solve for the derivatives, and then substitute a value for l, Solve won't choke on the algebra. The long time it takes is ...

Top 50 recent answers are included