# Tag Info

## New answers tagged numerics

7

Use Sinc[x] rather than Sin[x]/x BallooningFile = {0., 0.000136, 0.000572, 0.001152, 0.001907, 0.003004, 0.004199, 0.005479, 0.006834, 0.008256, 0.008985, 0.009738, 0.011271, 0.01285, 0.013651, 0.014468, 0.016119, 0.017797, 0.019496, 0.021211, 0.022069, 0.022934, 0.024661, 0.025522, 0.026386, 0.028103, 0.029806, ...

5

The PropertyValue::pvobj messages emitted are a bug, however they do not affect the functionality. The EquationTrekker GUI window should still open and operate normally (I've tried it on Windows, though Linux may have further problems). EquationTrekker is based on GUIKit which is being deprecated as far as I know, so this bug may not get a fix.

0

The following is a numerical solution to the toy model. It's too long to be a comment and not related to the question I really want to solve. Maybe someone can get help from the solution. $$\int_{-2}^{2}\int_{-2}^{2}\frac{1}{1-(x^2+y^2)}\, dx\, dy = 2\pi\int_{0}^{2}\frac{r}{1-r^2}\, dr+\int_{\mathrm{rest\ region}}\frac{1}{1-(x^2+y^2)}\, dV$$ In ...

3

When I first made this answer I was bleary eyed and didn't realize that Willinski's answer matched the question except that he made a nice edit by replacing, for example, 1.5 with 3/2. This answer is in addition to Willinski's fine work. I followed his procedure. I wanted to do a numerical study and try to find the region of interest. A = 1; J = 1; c = 1; ...

2

A = 1; J = 1; c = 1; beta = 1; int = Integrate[E^(c*beta*l*H*(3/2 x^2 - 1/2)), {x, 0, 1}] //FullSimplify; sum1 = Sum[l^(3/2)*E^(-u*l)*int, {l, 1, 100}]; sum2 = Sum[l^(5/2)*E^(-u*l)*int, {l, 1, 100}]; eq1 = A*beta^(-3/2)*E^(-J*beta - 1/2*c*beta*H^2); eq2 = eq1*(1 + sum1) - 1; eq3 = eq1*(1 + sum2) - H; ContourPlot[{eq2, eq3}, {u, 0, 10}, {H, 0, 5}, ...

3

Setting WorkingPrecision -> 5 in the gn integral gives you a reasonable convergence time. At the expense of some more computation time you can check that methods DifferentialEvolution and SimulatedAnnealing both return the sme result given here up to four decimal places. m1[x_] := 2 (x - 1)/x; m2[x_, y_] := x (2 - x y)/(2 (x - 1) y); fn[x_?NumericQ, ...

5

Your problem is partly that NIntegrate calls a function that at the given point in time is not full numerical (the first expression in fn). A way around that is to define the functions in such a way that hey will only evaluate for purely numerical values: Clear[m1, m2, fn, gn]; m1[x_?NumericQ] := (2 (x - 1))/x; m2[x_?NumericQ, y_?NumericQ] := (x (2 - x ...

6

Using polar coordinates r and f, the region of integration is given by { 0<= r <=2/Cos[f], 0<= f <= 2 \[Pi] } We can then proceed as follows. First integration with PrincipalValue g = 8 Integrate[r/(1 - r^2), {r, 0, 2/Cos[f]}, Assumptions -> 0 < f < \[Pi]/4, PrincipalValue -> True] (* Out[451]= -4 Log[-1 + 4 Sec[f]^2] *) ...

4

As mentioned in the comments and as the error message says, FindMaximum only accepts integer domain constraints for linear optimization problems, while of course nx * wx is a non-linear term. For ILP problems, FindMaximum uses a specialized solver from the COIN-OR branch-and-cut (CBC) library. NMaximize is using a different approach.

0

As already noted, due to the large range of variation between the nodes and weights of the Gauss-Laguerre rule, one would usually want to use arbitrary precision evaluation for high orders. In any case, let me present two alternative approaches to generating the nodes and weights for Gauss-Laguerre quadrature. I'll be linking to the papers explaining these ...

5

I think your problem is with the specification of the base case. You're basically trying to do the following, I think: f[n_] := f[n-1] + 1 f[5] /. f[0] -> 0 Of course, this fails because f[5] tries to evaluate completely before ReplaceAll even sees it. The solution is to specify a base case in the definition of the function: f[n_] := f[n-1] + 1 f[0] ...

5

Increase WorkingPrecision: NSolve[PDF[BinomialDistribution[80, p], 0] == 0.95200 && 0 < p < 1, p, Reals, WorkingPrecision -> 50] PDF[BinomialDistribution[80, p], 0] /. % (* {{p -> 0.00064096067673218860969986162632491931947341012861}} {0.9500000000000000000000000000000000000000000000000} *)

2

The use of NumericQ as mentioned by MarcoB and Guess who it is. in the comments seems to be important. Also, estimating the integral using a Total[Table[]] as in the code I posted in the second revision of the question makes this method computational feasible enough to solve my problem.

3

Following Szabolcs advice, it seems that the Documentation page for EvaluationMonitor contains all what you need for Method -> "Newton": data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}, {6, 4}, {7, 5}}; Clear[evalCount]; evalCount[_] = 0; nlm = NonlinearModelFit[data, Log[a + b x^2], {a, b}, x, EvaluationMonitor :> ++evalCount["Function"], Gradient ...

2

i = 2; j = 2; FixedPoint[Plus[#, f1[i++]/f2[j++]] &, f1[1.]/f2[1]]

4

You can find when $f_{2}(x) = 10^{-15}$, and then calculate the summation: f1[x_] := Exp[-x]; f2[x_] := 1/x; mybound = 10^-15; maxx = (x /. Solve[f2[x] == mybound, x])[[1]]; mysum = Sum[f1[x]/f2[x], {x, 0, maxx}] N[mysum, 10] (* 0.9206735942 *)

3

the NestWhile approach f1[x_] := Exp[-x]; f2[x_] := 1/x NestWhile[ {#[[1]] + f1[#[[2]]]/f2[#[[2]]] , #[[2]] + 1} & , {0, 1} , f2[#[[2]]] > 1/1000 & ] // First // N 0.920674 ( NSum is most certainly the better approach unless you have some peculiar functions )

7

Have you tried using NSum? In your question, it seems to me that you try to sum numeric values (instead of analytically evaluating a sum) and I think NSum is better for this. Simple example: f1[x_] := x; f2[x_] := Exp[x]; NSum[f1[x]/f2[x], {x, 1, Infinity}] (* 0.920674 *)

0

The analogous 1D problem, NDSolveValue[{D[u[x], x] == -u[x], D[v[x], x] == -v[x], u[0] == 1, v[1] == u[1]}, {u, v}, {x, 0, 2}] integrates across x == 1 without difficult. Therefore, I expected that NDSolveValue[{D[u[x, t], t] + D[u[x, t], x] == -u[x, t], D[v[x, t], t] + D[v[x, t], x] == -v[x, t], u[x, 0] == E^-x, v[x, 0] == E^(1 - x), ...

2

The error message is very descriptive: you haven't specified a boundary condition, but a condition on the middle of the region. If you amend the {x, 0, 2} to {x, 0, 1} it works correctly (but points out that your conditions are inconsistent, which they are). You should be able to impose conditions on the inside of the region by splitting the region in two, ...

3

The answer to the second part of your question is that you need to abandon machine precision arithmetic. Perhaps the best way to proceed is use Mathematica's exact arithmetic. a = Rationalize @ {15., 15.01, 3., 3.01} Round[Abs[a - Mean[a]], 1/100] {6, 6, 6, 6} You could also use Mathematica's slower but more accurate arbitrary precision arithmetic. ...

3

FullForm seems to be working fine here. Here is what belisarius alluded to in his comment, using your definitions and either InputForm, or FullForm as he suggested: InputForm@Abs[a - Mean[a]] FullForm@Abs[a - Mean[a]] (* Out from InputForm: {5.995000000000001, 6.005000000000001, 6.004999999999999, 5.994999999999999} Out from FullForm: ...

4

Using the formula given in the arXiv preprint Patrick linked to for the "carefree constant" gives: Exp[NSum[(-1)^k PrimeZetaP[k] (1 - LucasL[k])/k, {k, 2, ∞}, Compiled -> False, Method -> "AlternatingSigns", NSumTerms -> 20, WorkingPrecision -> 30]] 0.704442200999165592738713909247 Note that this agrees with the result in the OEIS ...

3

As it turns out, one can exploit the behavior of Interval[] when applied to a machine-precision number to obtain the previous and next representable machine-precision numbers (thanks to Szabolcs for the fix): SetAttributes[nextafter, Listable]; nextafter[x_?MachineNumberQ, s_?NumericQ] /; s != 0 := First[Interval[x]][[ -Sign[s - x] ]] To obtain a ...

0

You can use the finite element method with the method of lines as @toadatrix suggested, but for the FEM method to work, you need to do a little more. The Neumann boundary conditions need to be specified using NeumannValue. h[x_] := x*(30 - x)/900; op = D[u[t, x], t] - D[u[t, x], x, x]; begin = 0; end = 30; bc = {u[0, x] == 100*h[x]}; neumann = ...

2

So following @Szabolcs advice, {nl,ncounts}=Block[{c = 0}, {NonlinearModelFit[data, Log[a + b x^2], {a, b}, x, Method -> "Newton", EvaluationMonitor :> c++], c}] or (doing something slightly different) counting the number of Log Calls. Trace[nlm = NonlinearModelFit[data, Log[a + b x^2], {a, b}, x, Method-> ...

8

I think there's a bug in the internal function NDSolveSPRKDumpCheckSeparability that leads NDSolve to conclude that the system is not separable. I think you should report it and see if WRI can verify it (they would probably appreciate a link to this Q&A). It's a fair amount of work to track it down, and there is a lot of nearly unreadable stuff to ...

5

Erfc[-30. + 10^-1 I] used to return the result shown in the documentation through version 7.0.1. The implementation changed for version 8.0 and it started giving a machine precision answer (which is correct, more consistent and still demonstrates the same possible issue by being very close to 2). The (documentation) bug is that this example did not get ...

4

As Daniel Lichtblau showed in his comment, use exact numbers (or Rationalize) for input values alpha = 639/100; beta = 369/100; a1 = (1 - alpha)/2; a2 = (2 - alpha)/2; a3 = (1 - beta)/2; a4 = (2 - beta)/2; a5 = 1; b1 = 0; b2 = 1/2; SNR = 0; SNR0 = 10^(SNR/10); z = 2*(SNR0/(alpha*beta))^2; p1 = 2^(alpha + beta - 3)/(Pi*Sqrt[Pi]*Gamma[alpha]*Gamma[beta]); ...

4

I believe this has to do with the fact that intervals "grow" just a bit on evaluation with machine numbers to ensure that values at the endpoint will be included in the interval. NestList[Interval @@ # &, Interval[{0., 1.}], 5] // InputForm (* {Interval[{-2.2250738585072014*^-308, 1.0000000000000002}], Interval[{-4.450147717014403*^-308, ...

1

[Personally, I would be satisfied with SetAccuracy, at least in the use-cases in which I imagine I would need it. It just seems easier to me to learn how to work with the system instead of around it. Nonetheless, it seems to be possible....] Here's an idea of what I was talking about with \$PreRead in a comment. On a syntax error, it might fail ...

0

Perhaps you can use the Notation package to help. << Notation` Notation[a_ \[ScriptA] b_\[DoubleLongLeftRightArrow] N[Rationalize[a_,0],{Infinity,b_}]] I used a picture of a cell here to make it clear that this has been entered via the Notation palette Now you can use list=Table[10.75\[ScriptA] k,{k,1,10}] (* ...

4

Already answered in the comments by DumpsterDoofus and Daniel Lichtblau, to summarize: Machine floating point numbers such as 0.2 are not always exactly representable in binary (no terminating expansion in base 2). Thus floating point arithmetic is susceptible to roundoff error and other accuracy problems. For example, the following are not exactly equal to ...

Top 50 recent answers are included