# Tag Info

34

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

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

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

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

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 ...

10

I'll change your proposed function because it's a constant and for such a function all sums will be equal.So: f[x_] := 1/4 x^2 di = (2 - 1)/50; intervals = Range[1, 2, di]; leftSum = Sum[f[i] di, {i, Most@intervals}]; rightSum = Sum[f[i] di, {i, Rest@intervals}]; middleSum = Sum[f[i + di/2] di, {i, Most@intervals}]; exact = Integrate[f@x, {x, 1, 2}]; ...

10

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

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

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 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 Besides trivial observations that one cannot evaluate numerically integrals involving symbolic constants there are more interesting aspects of the problem at hand. First one should realize that a standard numeric approach is not appropriate for this kind of problems, since the integrand involves singular points (zero in the denominator) thus it is not ... 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 As of Version 9, you can work with vectors in NDSolve[]!: gravity = 10; withDrag[p0_, v0_, drag_] := Module[{p}, p[t_] := Evaluate@Array[Unique[][t] &, 3]; p[t] /. NDSolve[{ p[0] == p0, p'[0] == v0, p''[t] == drag*Norm[p'[t]]*p'[t] + {0, 0, -gravity}}, p[t], {t, 0, 5}] // ... 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}] == ...

6

Just consider the first integral. expr = x (A Ac (m2 + m1 (-1 + x)) + V Vc (m1 + m2 - m1 x))/(8 π^2 (-mh^2 (-1 + x) + (m2^2 + m1^2 (-1 + x)) x)); denominator = Collect[Denominator[expr], x] 8 mh^2 π^2 + 8 (-m1^2 + m2^2 - mh^2) π^2 x + 8 m1^2 π^2 x^2 It has two singular points. sol = Solve[denominator == 0, x] // Simplify If the singular ...

6

In Version 10, once the points have been obtained as per user21's approach, we can tetrahedralize them directly using DelaunayMesh pf = {Cos[u], Sin[u] + Cos[v], Sin[v]}; pp = ParametricPlot3D[pf, {u, 0, 2 Pi}, {v, -Pi, Pi}] data = Reap[ParametricPlot3D[Sow[pf], {u, 0, 2 Pi}, {v, -Pi, Pi}]][[2, 1]]; pts = Cases[data, {_?NumericQ, _?NumericQ, ...

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