9

We construct an NDSolve method which can pass an NIntegrate method to NIntegrate to set up an integration rule. We define a method nintegrate implements such a method. The requirements are the ODE is of the form y'[x] == f[x], and the NIntegrate method returns an interpolatory rule. Example: foo = NDSolveValue[{y'[x] == Sin[x^2], y[0] == 0}, y, {x, 0, 15}...


9

The line Table[N[PHI30EQ[u], 30], {u, 0, 0.1, 0.001}] doesn't do what you think it does. You're asking for 30 digits of precision, but you supply u as a machine number. If you mix arbitrary precision and machine precision like that, you'll get machine precision answers. I suspect you instead want is: Table[N[PHI30EQ[u], 30], {u, 0, Rationalize[0.1], ...


8

First of all, there do not seem to be any singularities in the integration region: $Assumptions = And[-Pi<=x<=Pi,Pi<=y<=Pi,Pi<=z<=Pi,-5<=r<=2]; den = 478 + (96*I)*br + 32*br^2 + 96*r - (64*I)*br*r - 32*r^2 + 11*Cos[2*x] - 264*Cos[y] - (48*I)*br*Cos[y] - 48*r*Cos[y] + 11*Cos[2*y] - 336*Cos[z] + 144*Cos[y]*Cos[z] + 12*Cos[...


6

We can integrate in 3 steps: Integrate[(yp/(b + yp^2)^(3/2)), {yp, 0, Infinity}, Assumptions -> b > 0]*(x - xp) /. {b -> g (x - xp)^2} //Simplify Out[]: (x - xp)/Sqrt[g (x - xp)^2] So we have intyp=1/Sqrt[g] as results and it means that Q[g,x] not depends on x. Next step: Integrate[(Exp[-ypp^2/(8 T)])*(-1 + ypp^2/(8 T)) (ypp/(g*(xp - ...


4

The first NIntegrate spends some time doing symbolic pre-processing of the integrand. You can turn that off and the integrals complete in about the same time: NIntegrate[ Exp[2 I s] Exp[2 I t] ((Cos[s] - Cos[t])^2 + (Sin[s] - Sin[t])^2) + 1, {s, 0, 2 π}, {t, 0, 2 π}, Method -> {Automatic, "SymbolicProcessing" -> 0}] // Timing (* result: {...


4

This interesting problem can be solved numerically by computing InterpolationFunctions for the two sums of integrals in the last two lines of code in the Question. λ2 = -(1/2) + r + 1/(1 + r) - r^2 Log[1 + 1/r]; k = 0.5; r = 0.5; t = Flatten[Table[ R1 = ImplicitRegion[θ (v + r) > r && 2 θ v s + (1 - θ) r > λ1, {{θ, 0, 1}, {v, 0, 1}}]; R2 = ...


4

I want to focus on the numerics here. There is a way to do this integral taking infinite number of singularities into account. One can split the integral into the domains containing only one singularity at $a x_n=\pi n$, i.e. $x\in[x_n-\frac{\pi}{2a},x_n+\frac{\pi}{2a}]$ and an integral in the interval $x\in[0,\frac{\pi}{2a}]$. On the last step we use NSum ...


3

Using the substitution u==a x the identity is transformed to Integrate[u/((u^2 + ab^2) Sin[u]), {u, 0, Infinity}] == Pi/(2 Sinh[ab]) The integral is singular at u==k Pi, k=1,2,... and the integration range is splitted into subintervals (similar to @yarchik answer) containing only one singularity (thanks to @flinty and @ ChipHurst comments): int[ab_?...


3

Following the @flinty's comment, we obtain Residue[x/(x^2 + b^2)/Sin[a*x], {x, π/a*n}, Assumptions -> n ∈ Integers] (*((-1)^-n n π)/(a^2 b^2 + n^2 π^2)*) Sum[%, {n, -∞, ∞}] (*0*) and 2*π*I*Residue[x/(x^2 + b^2)/Sin[a*x], {x, I*b}] (*π Csch[a b]*) Now, making use of the Jordan's lemma, we conclude that $$PV\int_{-\infty}^\infty\frac x {(x^2+b^2)\sin(...


3

Here is a method for estimating the error of integration. By sampling at two different rates, one can compare the integrals and estimate the error from the degree of precision of the integral (equals the InterpolationOrder in the cases at hand). If we interpolate $f(x)$ by a polynomial $p(x)$, then the error $\int_a^b f(x)\,dx - \int_a^b p(x) \, dx$ is ...


3

The symbolic processing misses the obvious way to compute these integrals, namely, Method -> "Trapezoidal", about 100 times faster than "SymbolicProcessing" -> 0: NIntegrate[ Exp[2 I s] Exp[2 I t] ((Cos[s] - Cos[t])^2 + (Sin[s] - Sin[t])^2) + 1, {s, 0, 2 π}, {t, 0, 2 π}, Method -> "Trapezoidal"] // RepeatedTiming ...


3

The integrand can be simplified to a significantly more compact form. First, the floats in f[x,y,z] can be replaced by rationals, which are easier for Mathematica to deal with: f[x_,y_,z_] = {{3/2 - (3*(Cos[x] + Cos[y]))/4, (-I)*Sin[x] - Sin[y], 7/2 - 3/(2*E^(I*z)) - (3*(Cos[x] + Cos[y]))/2, 0}, {I*Sin[x] - Sin[y], 3/2 - (3*(Cos[x] + Cos[y]))/4, 0, ...


3

J = ImplicitRegion[ Abs[xa] < 1 && Abs[xa - xb] < 1 && Abs[-L + xb] < 1 && Abs[xa - xc] < 1 && Abs[-L + xc] < 1 && L - 2 <= xa <= 1 && L - 1 <= xb <= 2 && L - 1 <= xc <= 2, {xa, xb, xc}]; Assuming[2 < L < 3, RegionMeasure[J]] (* 1/3 (27 - 27 L + 9 L^2 - L^...


3

reg = ImplicitRegion[{Abs[xa] < 1, Abs[xa - xb] < 1, Abs[-L + xb] < 1, Abs[xa - xc] < 1, Abs[-L + xc] < 1, L - 2 <= xa <= 1, L - 1 <= xb <= 2, L - 1 <= xc <= 2}, {xa, xb, xc}]; Assuming[2 < L < 3, Integrate[1, Element[{xa, xb, xc}, reg]]] // Simplify // AbsoluteTiming (*{0.142462, -(1/3) (-3 + L)^3}*)


3

The root cause of the difficulties encountered by the OP is that MeijerG, as implemented in Mathematica, often is not accurate when evaluated at machine precision. For instance, with the parameters, params = {p -> 1, q1 -> 1, q2 -> 1, q13 -> 1, qi1 -> 1, mu -> 3/2, L -> 3, gammak1 -> 1, BI -> 9, C1 -> 4, Bsd -> 4, M -&...


2

In this code we can check GaussianQuadratureWeights and FindRoot for potential errors. Let us evaluate GaussianQuadratureError[2 M, (1 + z) Sin[u[x/2 (z + 1)]], -1, 1] and we have answer for $u(x)=x^2$ -6.5402263142525195*^-105* Derivative[64][(1 + z)*Sin[(1/4)*x^2*(1 + z)^2]] Since $-1\le z\le 1, 0\le x\le 1$ we can conclude that Gauss quadrature not ...


2

I now realize that the question can be solved analytically, for the most part, although not with ImplicitRegion. The constraints embodied in R1 and R2 can be solved to obtain θ in terms of v and parameters. r1c1 = Reduce[θ (v + r) > r && v > 0, θ] // Last (* θ > 1/(1 + 2 v) *) r1c2 = Reduce[2 θ v s + (1 - θ) r > λ1 && v > 0 &...


2

Based on the detailed answer by @Hausdorff and the hints from @Akku14, I found that it is more efficient to transform the space integral into ParallelSum and then perform NIntegrate at each point with respect to r. Now the whole process takes only 15 min for all steps and with acceptable accuracy as shown in the attach Fig. Note that the max. value is about ...


2

Too long for a comment. First, your F[w_,s_]:= (1/(s + (I*a*w)/mx + ((\[Sigma]^2 + a^2)*w^2)/(2*mx) + balpha/mx^2*s^(alpha - 1)*(I*w*a + 1/2*(\[Sigma]^2 + a^2)*w^2))) is a complex valued function. I don't see any reason why its inverse Fourier transform should be a real valued function. Second, up to the definition used by Mathematica, its inverse Fourier ...


2

Updated to Accommodate Change in $a$ Definition The OP changed the definition of the symbol, $a$ from my original answer. The code is essentially the same as before: tf = 10; dt = 0.1; N1 = 10; Ω = 2; a = ConstantArray[0, N1 - 1]; a[[N1 - 1]] = 1; zeros = ConstantArray[0, N1 - 1]; Y1[t_] := Table[Subscript[x, i][t], {i, 2, N1}]; F[t_] := Table[Subscript[f, ...


2

We can use first order approximation on z as follows d = 84.9601; tau = 0.1; t0 = .1; OmegaC = 27.7259; pulseTime = 2 Pi/OmegaC; storageTime = 0.2; twriteon = t0 - pulseTime/2; twriteoff = t0 + pulseTime/2; treadon = t0 + pulseTime/2 + storageTime; treadoff = t0 + pulseTime/2 + storageTime + pulseTime; Clear[OmegaFunc] OmegaFunc[t_] := Piecewise[{{0, t &...


2

My guess is NIntegrate does not recognize the singularities of HeavysideTheta. At least, if you specify it, you get rid of the problem. NIntegrate[y x^2 HeavisideTheta[5 - x - y], {x, 0, 10}, {y, 0, 10}, Exclusions -> 5 - x - y == 0, AccuracyGoal -> 16] (* 52.0833 *) The default setting for AccuracyGoal is Infinity, which means NIntegrate will ...


2

Try this: NIntegrate[Exp[-(x^2 + y^2) - z]/ Sqrt[x^2 + y^2], {z, 0, 10}, {x, -5, 5}, {y, -5, 5}] (* 5.56808 *) Have fun!


2

You had some syntax errors from your brackets which should be parentheses and also with Log( ... ) which should be Log[ ... ]. The following works: h = 1 b = Pi/(2 h) g[w_?NumericQ] := -(2 Sinh[b] Cosh[b])/(h (Sqrt[Cosh[Pi w/h] - 1] (Sqrt[Cosh[Pi/h] - Cosh[Pi w/h]]))) r1[w_, x_, y_] := Log[Cosh[(w-x) Pi/(2 h)] + Sin[(h-y) Pi/(2 h)]]/(Cosh[(w-x) Pi/(2 h)] ...


1

This should do it: var = Array[x,n] sol = NIntegrate[f@@var, Evaluate[Sequence @@ ({#, -Infinity, Infinity} & /@ var)]] Test it: n = 5; f = (Exp[-#1] + Exp[-#2] + Exp[-#3] + Exp[-#4] + Exp[-#5]) & var = Array[x,n] sol = NIntegrate[f @@ var,Evaluate[Sequence @@ ({#, -1, 1} & /@ var)]]


1

I rewrote this without the subscripts, I added NumericQ in places where it was needed, I used SetDelayed (:=) for some of your functions, and I assumed by i you meant I the imaginary number: Z = 6; n = 6; K1 = 1.55; fpn1[q_?NumericQ] := (18.25974896615874*(1 - 0.8748275119106319*I)*I)/E^(0.2115*q^2) Fpn1[q_?NumericQ] := ((4*Pi)/q)*NIntegrate[ρn[r]*Sin[q*r]*...


1

The integration NIntegrate[g[NIntegrate[f[x],{x,0,y}]],{y,0,1}]fails because the inner NIntegrate doesn't know about y to be numeric! Try int[y_?NumericQ]:=NIntegrate[f[x],{x,0,y}] NIntegrate[int[y],{y,0,1}]


1

If you you only need a numerical result for given L try int[L_?NumericQ] :=NIntegrate[ Boole[Abs[xa] < 1] Boole[Abs[xa - xb] < 1] Boole[ Abs[-L + xb] < 1] Boole[Abs[xa - xc] < 1] Boole[ Abs[-L + xc] < 1], {xa, L - 2, 1}, {xb, L - 1, 2}, {xc, L - 1,2} ] int[2.5] // AbsoluteTiming (*{0.0222419, 0.0416667}*)


1

This is not a complete answer but an extended comment. ImplicitRegion does not like the usage of function Tmaxc12, so we can construct it inline: region[U_?NumericQ] := ImplicitRegion[ (0.5*MC12 (vx^2. + vy^2. + vz^2.) + (1 - Cos[theta])*(Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2] + MC12*vz)* Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]/MC12 - Sqrt[...


1

On such a high dimensional integral, the default rule is the Monte Carlo rule. You can increase the number of points. I also increased the PrecisionGoal, so that the error estimate will be reported. NIntegrate[..., Method -> {"MonteCarloRule", "Points" -> 10^6}, PrecisionGoal -> 6] NIntegrate::maxp: The integral failed to ...


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