Hot answers tagged

24

The simplest way to make new NIntegrate algorithms is by user defined integration rules. Below are given examples using a simple rule (the Simpson rule) and how NIntegrate's framework can utilize the new rule implementations with its algorithms. (Adaptive, symbolic processing, and singularity handling algorithms are seamlessly applied.) Basic 1D rule ...


19

Motivation (for a new semi-symbolic integration strategy) Consider the following integral, which cannot be done neigther by Integrate: Integrate[BesselJ[y, x^3], {x, 0, ∞}, {y, 0, 1}] (* Integrate[If[Re[y] > -(1/3), Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]), Integrate[BesselJ[y, x^3], {x, 0, Infinity}, Assumptions -> Re[y] <= -(1/3)]...


7

As mentioned in the comments, here's the answer to how I got the median. I found the value $y$ that minimized $$ f(y) = \int_0^1 \left| 2^x-x-y \right| dx. $$ First, to break up the absolute value in the integrand, I solved for when $2^x = y$: Solve[(2^x - x) == y && 9/10 < y < 1, x, Reals] Then I integrated and found where $f'(y) = 0$: ...


6

It's not a bug and it's not so uncommon. For an explanation have a look here. This and some related issues also appear in this MathGroup thread. Also relevant: 1 2.


6

This is not an answer. But I don't believe we should close this question as "easily found in the documentation". Numerics in Mathematica is an extremely complicated and mostly undocumented subject, where several mathematical concepts run up against each other in subtle and non-trivial ways. I have been thinking for some time that we ought to address this ...


5

What is observed for testfunc and integrand2 is explained with the use of adaptive sampling, symbolic processing, and singularity handlers by NIntegrate. The home-cooked Simpson integration strategy in the question is too simple for these integrands. For the function integrand it seems that only the adaptive sampling gives the advantage of NIntegrate. Needs[...


4

Here is an expanded version of JM's implementation of your method that he proposed in comments. As an aside, I found this interesting write-up by Michele Benzi (Emory University) on Cimmino, his method and other accomplishments, and the Italian school of numerical analysis in the 1920s-30s; a very interesting read. I first propose a version using the ...


4

Here is another type of integration rule. It has a customized error estimator, incorporating which was my main interest. In fact, the rule below computes an integral of arbitrary dimension, using its own internal formula for the integral as well. Introduction A rule of the type NIintegrate`GeneralRule[{abscissas, weights, errweights}] computes the ...


3

The reason is that the Element function does not test the Wolfram representation of the number (exact integer, floating point real, etc.), but the mathematical meaning of the number (integer, real, etc.). Compare the followings: Element[1., Integers] (* 1. ∈ Integers *) Element[1 + 0. I, Reals] (* 1. + 0. I ∈ Reals *) These examples indicate that ...


3

They are not identical computations. With the first form, (mu/2 gt).gt Mathematica can take advantage of vector arithmetic, usually going through specialized routines like LAPACK. The second form, Sum[(mu[[i]]/2 gt[[i]]) gt[[i]], {i, Length@mu}] however, will usually be calculated term by term because there is a possibility that the input can change ...


3

Thanks to the recent post by ilian who uncovered and explained the undocumented Optimization`FindFit`ObjectiveFunction, I have found Optimization`FindFit`ResidualFunction which seems to be ideally suited to the original goals of the question. With the latter we don't even need to construct the residual vector explicitly, we just specify the model function ...


3

Once you have the NDEigensystems code in place, it is not difficult to extend it to work on multiple regions (this is, I believe, what @JM was nudging you towards). For instance, generate a list of such polygons using Table: r = Table[ Graphics[ Polygon[{ {0, 0}, {0, 1}, {1, 1}, {1, (1 + epsilon)/2}, {2, (1 + epsilon)/2}, {2, 1}, {3, 1}...


2

If you have a C compiler you could use C's nextafter directly. Needs["CCompilerDriver`"] ClearAll[nextafter] " #include \"WolframLibrary.h\" DLLEXPORT mint WolframLibrary_getVersion() { return WolframLibraryVersion; } DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) { return LIBRARY_NO_ERROR; } DLLEXPORT void ...


2

This can be done with "EquationSimplification" -> "Residual", provided the denominators are cleared from the equation. ndsol = NDSolve[{x y''[x] + y'[x] + x y[x] == 0, y[0] == 1, y'[0] == 0}, y[x], {x, 0, 10}, Method -> {"EquationSimplification" -> "Residual"}]; Check with DSolve: dsol = DSolve[{x y''[x] + y'[x] + x y[x] == 0, y[0] == 1, y'[...


2

I think that the imaginary component is a rounding error. Try N[exp, 40] and the imaginary part is returned as zero. ADDED If evaluated to sufficient significant figures, we see that the expressions given above all have an imaginary part indistinguishable from zero. N[exp // Simplify, 40] N[exp // ExpandAll, 40] N[exp // ExpToTrig, 100] N[exp // ...


2

I was looking for a previous question of which this might be a duplicate; although many previous discussions hinge on precision issues, and the difference has been mentioned in comments, I could not dredge up an explicit discussion of $MachinePrecision vs MachinePrecision through the SE search engine, so I thought I might sum up the discussion in comments, ...


2

In general one should not expect to obtain a general symbolic solution (a function x[a,b]) to the given equation since there are two independent variables, see e.g. Solve symbolically a transcendental trigonometric equation and plot its solutions for certain aspects regarding transcendental equations. To get an idea how the solution depends on parameters a ...


2

Perhaps this! It works if you only need the plot and not the values, or the function. That would require more work. ContourPlot3D[ Evaluate[-x + (Log[f2[x]] - Log[f1[x]] + Log[b] - Log[a])/((1 - a) + (1 - 1/f1[x]) - (1 - b) + (1 - 1/f2[x]))] , {a, 0, 1}, {b, 0, 1}, {x, 0, 1} , AxesLabel -> {"a", "b", "x"} , Contours -> {0} , Mesh -> None ...


2

Seems to be a case of ill conditioning of the input system. If I redo using exact input and set NSolve to work on high precision then I get a plausible outcome. c1 = 50; c2 = 3/2; c3 = 2/5; c4 = 1/2; c5 = 8; c6 = 16/5; k1 = 37*10^(-17); k2 = 83*10^(-8); k3 = 87*10^(-12); k4 = 32*10^(-6); k5 = 86*10^(-32); eqs = {f1*f2 - k1*x1^2, f1*x3 - k2*x2*x1, f1*(c6 ...


1

What happens here is that Association has attribute HoldAllComplete and hence does not evaluate x[1] to a numerical value: Clear[x]; assoc = <|"x" -> x[1]|>; x[1] = 0; assoc <|"x" -> x[1]|> Numerical functions work by assigning values to the variables and then evaluating the objective function, they don't perform replacement like in ...


1

This reminds me of NSolve finds real-valued results in version 9, but not in version 10, in which you can use any of the Method settings from Methods for NSolve "EndomorphismMatrix" "CompanionMatrix" "Legacy" "Aberth" "JenkinsTraub" or even a nonexistent method "Foo": NSolve[eqs, {x1, x2, x3, x4}, Reals, Method -> "Foo"] (* {{x1 -> 23.2915, x2 -&...


1

While not being a Mathematica expert, I assume that everything works as expected: your 1.1 "actually is" 1.1000000000000000888 since it is probably stored with a binary form in the IEEE754 standard (count the number of exact decimal digits and remember that IEEE754 standard with double-precision type has 15/16 exact decimal digits). See http://mathworld....


1

I can solve the toy problem with NMaximize[{x, Im[Sqrt[2 - x]] == 0}, x]


1

Try the following shory code: DeleteDuplicatesBy[lst,Floor[#,10^-4]&] Will this help?


1

It seems that the most straightforward way goes through RealDigits: nextAfter[x_Real] := FromDigits[MapAt[1 + # &, MapAt[0*# &, RealDigits[x, 2], 1], {1, -1}], 2] + x



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