Tag Info

Hot answers tagged

36

Control the Precision and Accuracy of Numerical Results This is an excellent question. Of course everyone could claim highest accuracy for her product. To deal with this situation there exist benchmarks to test for accuracy. One such benchmark is from NIST. This specific benchmark deals with the accuracy of statistical software for instance. The NIST ...


22

Solving 1D and 2D complex Schroedinger wave equations with NDSolve I do not agree with you when you write: I know the NDSolve is not magic... My opinion is that NDSolve is one of the most complex functionality I've met so far in the Mathematica environment, with its millions of options and special function this is a real complex thing and it is hard ...


15

According to the Mathematica documentation on this page: Here is how to define a 5(4) pair of Dormand and Prince coefficients [DP80]. This is currently the method used by ode45 in MATLAB. DOPRIamat = { {1/5}, {3/40, 9/40}, {44/45, -56/15, 32/9}, {19372/6561, -25360/2187, 64448/6561, -212/729}, {9017/3168, -355/33, 46732/5247, 49/176, ...


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


14

Some frames from my version of the animation: Here's the code I used: orbit[posStart_?VectorQ, derStart_?VectorQ] := Block[{c = -Rationalize[6.672*^-11*7*^17], x, y, z, t}, {x, y, z} /. First @ NDSolve[ Join[Thread[{x''[t], y''[t], z''[t]} == c {x[t], y[t], z[t]}/Norm[{x[t], y[t], z[t]}]^3], ...


14

I think it's worth pointing out that the problem can be solved "straightforwardly" (i.e., really using only NDSolve) once you know the options that Stefan used in ProcessEquations (which I upvoted because those options are the main ingredient): Below I show the original problem of a Gaussian wave packet with no initial momentum, and then a modified case ...


12

You can do the analytic integral in Mathematica too, by telling it to perform the upper integration limit as follows: With[ { i = Integrate[((r^3 - 7)^(2/3)*(1 - (r^3 - 7)^(2/3)/r^2))/r^3, r] }, Simplify[ Limit[i, r -> Infinity] - i /. r -> 2 ] ] (* ==> 23/64 + Pi/(3 Sqrt[3]) - 2 (-(1/7))^(1/3) Hypergeometric2F1[1/3, 1/3, 4/3, ...


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

You can express your integral in terms of a differential equation and use NDSolve. Since NDSolve builds up the solution as it goes, this is typically much faster. Clear[y]; y[x_] = y[x] /. First[ NDSolve[{y'[x] == Sin[x], y[0] == 0}, y[x], {x, 0, 10}] ]; t = AbsoluteTime[]; Plot[y[t], {t, 0, 10}] AbsoluteTime[] - t


11

According to the error message: NDSolve::nlnum1: "The function value {0} is not a list of numbers with dimensions {25} when the arguments are {50.,{<<25>>}." I think you should feed a 25-length list of 0 to n[x,t] in the WhenEvent: WhenEvent[t > 50, n[x, t] -> ConstantArray[0, 25]] Plot the result: Plot3D[Evaluate[n[x, t] /. sol], ...


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


10

Believe the numerical one. Mathematica simply could not do the symbolic integration. Symbolic integration will travel via a different code path. Here the symbolic integration done using Maple, and it agrees with the numerical solution given by Mathematica's NIntegrate The analytical answer is (7/18)*hypergeom([-1/3, 1, 1], [2, 2], ...


10

You can use Interpolation to construct a function from your data that can be passed to NIntegrate. Here's how: data={ {0.,0.,0.,0.,0.,0.,0.,0.}, {0.,-10.9421,-17.3061,-19.0045,-19.0045,-17.3061,-10.9421,0.}, {0.,-16.8109,-27.6012,-24.8577,-24.8577,-27.6012,-16.8109,0.}, {0.,-19.9862,-34.0245,-28.3369,-28.3369,-34.0245,-19.9862,0.}, ...


10

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


10

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


10

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


9

This integral equation is solvable using the LaplaceTransform technique: Clear[s, t]; eqn = y'[t] == -Integrate[y[t1] Exp[t1 - t], {t1, 0, t}] LaplaceTransform[eqn, t, s] (* ==> s LaplaceTransform[y[t], t, s] - y[0] == -( LaplaceTransform[y[t], t, s]/(1 + s)) *) Solve[%, LaplaceTransform[y[t], t, s]] (* ==> {{LaplaceTransform[y[t], t, s] -> ...


9

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


9

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


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


7

This is a classical shock-tube problem in which a initially diaphragm separates a hi-pressure, high-density region from one of lower pressure and density. The classical exact solution has multiple discontinuities, a shock wave and a contact-surface (density discontinuity) that propagate to the right, and a continuous rarefaction wave traveling into the ...


7

This question appears to be about Options, not optional arguments. While it is true that OptionsPattern[] is a form of optional argument that is not typically how options themselves are described. Edit: I see that the tutorial refers to options as "named optional arguments" so I guess my understanding of convention was incorrect. The meat of your question ...


6

An alternative approach is to form an approximate interpolating function from your actual function, which is (often) cheaper to integrate. To wit, nsin = FunctionInterpolation[Sin[x], {x, 0, 10}]; ni = Derivative[-1][nsin]; Plot[ni[t], {t, 0, 10}] In this case, we know that the integral of sine can be expressed simply, so we can compare the exact ...


6

As has been noted by ruebenko in the comments, there does seem to be a bug in the handling of infinite-range Bessel function integrals when MinRecursion and MaxRecursion are both set to non-default values. For instance, even the simple NIntegrate[BesselJ[0, x], {x, 0, ∞}, MinRecursion -> 10, MaxRecursion -> 15] chokes with a NIntegrate::minmax ...


6

I think your initial condition is singular. In order to solve the ODEs, consider $(x(t), y(t))$ as a planar curve, we may try changing the parameter $t$ to the arc length parameter $s$: $$\left\{\begin{split} \frac{\mathrm{d}x}{\mathrm{d}s}=\frac{x'(t)}{\sqrt{x'(t)^2+y'(t)^2}}\\ \frac{\mathrm{d}y}{\mathrm{d}s}=\frac{y'(t)}{\sqrt{x'(t)^2+y'(t)^2}}\\ ...


6

Since no comment or answer has given as of the time I saw this post, and I don't have enough reputations to leave a comment, let me give a quick answer I used to solve the same kind of problems. I'm sure there must be a better way to do it with Mathematica, so this is just a beginning. To evaluate an N-dimensional integral with a highly oscillatory ...



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