# Tag Info

29

The only reason I am attempting to answer this is to perhaps get a Reversal badge. There you go... We will go slowly and this answer is the basis for what comes next. Let's start with two dimensions. You'll see why. We create a rectangular region: Needs["NDSolveFEM"] mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}, "MeshOrder" -> 1, ...

16

After a lengthy study (I'm using version 8) I conclude that there is a bug in Mathematica in the Integrate function when applied to a Sqrt integrand. Ok. let's go (some patience is required because of the long text) Let us define the functions corresponding to your integrals. Remark: because of the relation $1 + cos(2x) = 2 cos^2(x)$ the two forms of ...

15

Here's my attempt. To get the matrix representing the Laplacian I use LaplacianFilter on an array of symbols and CoefficientArrays to extract the coefficients. n = 200; shape = ArrayPad[ConstantArray[0, {n/2, n/2}], {{0, n/2}, {0, n/2}}, 1]; shapeVector = Flatten @ Position[Flatten @ shape, 1]; symbolArray = Array[x, {n, n}]; symbolLaplacian = ...

12

In this second answer I give the cause for the mismatch in the integrals, show how to remove it, and make a suggestion to improve the function Integrate[]. Simplified restatement of the problem In order to focus on the core of the problem we consider the simpler integral $\int_0^1 \sqrt{\cos (2 π k r)+1} \, dr$. It has the square root and the cosine ...

12

While the other answers are nice, the icon deserves a closer look: Note, in particular, that four of the six edges are not constrained by the ostensible Dirichlet boundary conditions, nor is it clear that they solve a Neumann problem. And indeed, as I noted in the comments this is supported by the OP's first link. In short, to produce the logo, they took ...

11

This does not completely answer the question, but you can get some useful information from the undocumented option IntegrationMonitor. For example: NIntegrate[Sin[Sqrt[x]], {x, 0, 1}, IntegrationMonitor -> Print] You can see (in the Experimental`NumericalFunction) that the change of variables $\sqrt{x}\to x$ has been used to convert the integrand to ...

11

Why the order could matter Superficially, picking some values for t, integrating τ, and finally integrating the results over t is a different calculation than calculating it in a different order. Now there are a few things to explain and investigate to show that this naive observation has a bearing on the OP's integral. Broadly, I would say that the ...

11

The answer is no because of fundamental mathematical limitations which origin in the set theory regarding countability (see e.g. Cantor's theorem) - functions over a given set are more numerous than its (power) cardinality. Neither Mathematica nor any other system can integrate every function in even much more restricted class, namely Riemann integrable ...

11

I had this laying around from a course in numerical linear algebra I taught a few years ago. Here's a matrix whose nonzero elements describe the basic shape. size = 50; nw = Partition[Table[i, {i, 1, size^2}], size]; sw = Partition[Table[i, {i, size^2 + 1, 2*size^2}], size]; se = Partition[Table[i, {i, 2*size^2 + 1, 3*size^2}], size]; L = ...

11

You can write your own algorithm and use it from NDSolve. For example, for RK4: CRK4[]["Step"[rhs_, t_, h_, y_, yp_]] := Module[{k0, k1, k2, k3 }, k0 = h yp; k1 = h rhs[t + h/2, y + k0/2]; k2 = h rhs[t + h/2, y + k1/2]; k3 = h rhs[t + h/2, y + k2]; {h, (k0 + 2 k1 + 2 k2 + k3)/6}] CRK4[___]["DifferenceOrder"] := 4 CRK4[___]["StepMode"] := Fixed ...

10

Try not to supply machine numbers to integrals over infinite domains. They can cause errors that build up to the extent you have seen. Either compute the symbolic integral with exact numbers (and then convert it to a numeric value) L = 2; sn = 1; a = 10^(sn/10); b = 10^(sn/10); c = a/100; result = 2*Sqrt[1/Pi]*Integrate[(1/(E^z*Sqrt[z]))*(1 - (a/(a + ...

10

Just to contribute to the debate, here is some more evidence that supports the proposition that numerical error is the issue. If we run the integral through various permutations of the ways of making exact and approximate calculations, the pattern I think suggests that numerical error is the reason the OP's integral is so far off. (* the integrand and ...

10

Update Almost ten times faster again, or about 90 times faster than the OP's way (0.069 sec v. 5.46 sec): For the second integral, we can find its derivative with respect to x and then integrate with NDSolve. The derivative of the integral has two components, one from differentiating under the integral sign dxdz1 and one from plugging in the limit of ...

9

Update The problem is subtler than my first analysis revealed. There is indeed a problem with the variable et in NIntegrate not being properly blocked. Part of the problem has to do with the extra braces in firstFuncK which has the form {{f -> InterpolatingFunction[<>]}} Somehow that leads to an evaluation of et in the integrand f[et, k] /. ...

9

In Mathematica 10, this computation may be made as follows: Clear @ r volSphere9[r_] = RegionMeasure[Ball[ConstantArray[0, 10], r]] (π^5 r^10)/120 volSphere9[1000.] 2.55016*10^30

9

Use NDSolve antiD = NDSolveValue[{f'[x] == Sqrt[1 + x^3], f[0] == 0}, f, {x, 0, 10}] Example usage: Plot[antiD[x], {x, 0, 10}] Alternatively... This works because this function can be antidifferentiated (by Mathematica). antiD = FunctionInterpolation[ Evaluate @ Integrate[Sqrt[1 + x^3], {x, 0, t}, Assumptions -> 0 < t < 10], {t, 0, ...

9

The interpolation step seems to be unnecessary because the integral into which it enters can be equally well approximated as a Riemann sum. So to get really fast results you could do the following: {r, h} = Transpose[hrdata]; d = Differences[r]; Clear[s]; s[q_] := (4 Pi 0.83 )/q Total[d Rest[r Sin[q r] (h - 1)]] ans2 = Table[{i, s[i]}, {i, 0.05, 11.4, ...

8

Every integral over a function behaving asymptotically (when $x$ goes to infinity) as $\frac{1}{x^\alpha}$ where $\alpha \leq1$ is divergent, it's a mathematical theorem which could be found in every reasonable handbook of calculus. Since Tanh[ π Sqrt[x]] goes to one rapidly we find that the integral is indeed divergent. We can demonstrate this fact with ...

8

It's numeric integration. So it has no means of "knowing" the correct result is zero. In the process error estimates will be formed and if they are larger than the estimated result, this is a problem. But of course they must be larger since the actual result is zero. The way to tame this is to specify an AccuracyGoal that is attainable using the given ...

8

There is a useful attribute, NHoldFirst whose purpose is to protect the function from exactly that. So setting: SetAttributes[a, NHoldFirst]; and then evaluating the integral works the way you want: Integrate[(a[1] + x)^2, {x, 1., 2.}] (*2.33333 + 3. a[1] + 1. a[1]^2*) The relevant example from the documentation cites "indexed" functions that are ...

7

Here is a way to get an approximate symbolic expression for the integral. Some of the coefficients are approximate because at some points in the process we need the value of a definite integral at a = 1, and unfortunately those integrals have to be computed numerically. First, a lemma: Let $f(a) = \int_1^\infty g(x, a) \; dx$. Then $f(a) = f(1) + ... 7 Here is another way that uses the Graphics object directly: gr = ParametricPlot3D[{Cos[u], Sin[u] + Cos[v], Sin[v]}, {u, 0, 2 Pi}, {v, -Pi, Pi}] We discretize the graphics using DiscretizeGraphics mr = DiscretizeGraphics[Normal[gr /. (Lighting -> _) :> Lighting -> Automatic]] We compute the convex hull hull = ... 7 The sum of the squares should be less than or equal to r^2 rather than r. d = 10; r = 1000; F = Piecewise[{{1, Sum[x[i]^2, {i, d}] <= r^2}}, 0]; NIntegrate[F, {x[1], -1000, 1000}, {x[2], -1000, 1000}, {x[3], -1000, 1000}, {x[4], -1000, 1000}, {x[5], -1000, 1000}, {x[6], -1000, 1000}, {x[7], -1000, 1000}, {x[8], -1000, 1000}, {x[9], -1000, 1000}, ... 7 Here is the same approach as Jens, but using Association and assuming uniform grid data: st[hr_?MatrixQ, rho_, min_, max_, step_] := Module[{h = <|Rule @@@ hr|>, s, freq, dr = hr[[2, 1]] - hr[[1, 1]]}, s = Function[q, 4 Pi rho /q Tr[dr # Sin[q #] (h[#] - 1) & /@ Keys[h]]]; freq = Range[min, max, step]; Transpose[{freq, s /@ freq}] ] Now: ... 6 There are two ways come to my mind to go. 1. Truncate the upper limits: Since OP has an exponential decay term like what it reads, truncating the integral limits at 2000000 should give a reasonably precise result: NIntegrate[( E^(-(x - 1000000)^2/(2*200000^2)) *E^(-(y - 1000000)^2/(2*200000^2)) *E^(-(z ... 6 There is a way to get Mathematica to calculate the equation of a transformed cylinder, which can then be used to calculate the volume. First, since you're translating the cylinder, too, I rewrote your transformation to include the translation. We can also define inequalities to define the cylinder. xform[x_, y_, z_, a_, b_] := RotationTransform[a Pi/2, ... 6 The analytic answer is $$%\sum_{n=1}^\infty \frac{1}{2n(n+1)^2}\sum_{n=1}^\infty \frac{1+3n}{6n^2(n+1)^3}+\sum_{n=1}^\infty \frac{3n^2-1}{6n^2(n+1)^3}=\\ 1+\pi ^2\frac{2 \zeta (3)-9}{72} \approx 0.0958502$$ Therefore, Mathematica is correct. Proof The 3D integral NIntegrate[FractionalPart[x/y] FractionalPart[y/z] FractionalPart[z/x], {x, 0, 1}, ... 6 Using the functions defined in my answer to your previous question here you have the intersections with all three coordinate planes: getOneCluster[pts_List, maxDist_?NumericQ] :=(*Returns a cluster*) Module[{f}, f = Nearest[pts]; FixedPoint[Union@Flatten[f[#, {Infinity, maxDist}] & /@ #, 1] &, {First@pts}]] clusters[data_] := Module[{f, ... 6 Here's a way to get an exact symbolic expression for the integral. Use the notation fi[m,n,a] to denote Integrate[(Exp[-a x](x^2 - 1)^(m/2))/x^n, {x, 1, Infinity}]. Define integration by parts, where the boundary term vanishes. fi[m_, n_, a_] /; m > 0 && n > 1 := (m fi[-2 + m, -2 + n, a] - a fi[m, -1 + n, a])/(-1 + n); See how we're ... 6 This is already bult into Mathematica$Assumptions = {s > 0}; dist = NormalDistribution[m, s]; pdf[x_] = PDF[dist, x]; cf[w_] = CharacteristicFunction[dist, w]; pdf[x] == InverseFourierTransform[cf[w], w, x, FourierParameters -> {1, 1}] // Simplify True cf[w] == FourierTransform[pdf[x], x, w, FourierParameters -> {1, 1}] == ...

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