# Tag Info

8

The following is not particularly fast and could be more accurate but does make progress toward the goals set in the Question. To begin, consider the analytical properties of f, as defined in the Question. By observation, it has branch points at {I Sqrt[6/5] ξ, I Sqrt[2/5] ξ} and their conjugates. Poles are obtained by p /. ...

0

T[x_] := Piecewise[{{1 - x, 0 <= x < 1/7}, {(x + 6)/7, 1/7 <= x <= 1}}]; x[n_] := (1 - a[n - 1]) x[n - 1] + a[n - 1]*T[(1 - b[n - 1]) x[n - 1] + b[n - 1] T[x[n - 1]]]; x[0] = 0.1; a[n_] := n/(n + 1); b[n_] := n/(n + 5);

2

I'd suggest a simpler approach that seems to work just fine: Use arbitrary precision input if you can. In your case, instead of d -> 0.001/5, use d -> 5/1000 integrand = Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c - 1/2 a^2) (Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) /. {b -> 5, ...

3

The problem is that you are not using large enough precision goal for an oscillatory integrand with very small absolute values over the integration range. Look at the log plot -- there are more oscillations than the ones visible with Plot: The remedy is to use higher precision goal. Try this: NIntegrate[ Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/ 2 ...

2

Short answer: One workaround is to use Method -> "LevinRule" instead. Long answer: As mentioned by Xavier in a comment above, changing BesselJ[0, x] to BesselJ[0, Re[x]] resolves the issue: NIntegrate[{1, -1} BesselJ[0, Re@x], {x, 0, ∞}, WorkingPrecision -> 32, Method -> "ExtrapolatingOscillatory"] Precision /@ % ...

2

Is there any reason you can't solve for logarithms? If you do, NSolve handles the provided example pretty straightforwardly: Thread[y -> (Exp[y] /. NSolve[-y == Log[10^(-50)], y])] (* {y -> 1.*10^50} *)

6

The problem seems to be that you give the initial condition at a different point that your lower limit of the integration. If you change x1 by 0 things seem to work: (*INPUT*)MP = 1; m = 0.2 MP; c = Sqrt[3/2]; x1 = 0; x2 = 4000; (*ODE*) Clear[s2] Table[s2[i] = NDSolve[{y''[x] + c (Sqrt[m^2 y[x]^2 + y'[x]^2]) y'[x] + m^2 y[x] == 0, y[0] == ...

3

If all your problems are of this sort, simply use Solve. Solve acts by manipulating the symbols, which in these cases is trivial, so that Solve[1/x == 10^(-50), x] gives the proper answer directly. Mathematica recognizes this equation as $1/x = a$ and solves, yielding $a = 1/x = 10^{(-(-50))} = 10^{50}$. In contrast, NSolve performs a numerical ...

3

Analysis of the error (bug?) We can see from the trace below that the second limit, which carries out a ratio test for the product, mistakenly yields -17 (which would indicate divergence, if correct). Trace[ NProduct[(n^2)!/stirling[n^2], {n, 1, Infinity}], _Limit, TraceInternal -> True, TraceForward -> True] There might have been some ...

2

A point has measure 0 in any numeric approximation possible, so it's not only mathematica who doesn't like point-like conditions. soln = NDSolveValue[{Laplacian[u[x, y], {x, y}] == NeumannValue[-1, y == -1. && Abs[x] < 0.1] + NeumannValue[1, y == 1 && Abs[x] < 0.1], u[0, y] == 0}, u, {x, y} \[Element] Rectangle[{-1., ...

2

In version 10.2 or later you can use the EXPERIMENTAL function FindFormula lst = {{0, 1}, {5/999, 0.999925}, {5/333, 0.9997}, {25/999, 0.999325}, {35/999, 0.998801}, {5/111, 0.998127}, {55/999, 0.997305}, {65/999, 0.996335}, {25/333, 0.995217}, {85/999, 0.993952}, {95/999, 0.992542}, {35/333, 0.990986}, {115/999, 0.989287}, {125/999, ...

4

First of all there is somewhere in Mma a package for numerical calculation of derivatives, but I did not manage to find a reference. To offer a way to calculate the derivative. You could use the interpolation function. Here is your list: lst = {{0, 1}, {5/999, 0.999925}, {5/333, 0.9997}, {25/999, 0.999325}, {35/999, 0.998801}, {5/111, 0.998127}, ...

1

This is just a long comment really. Out of curiosity, I compared the result of three packages: Mathematica, MATLAB and eig_sym from Armadillo (compiled on OS X). You said that your C++ code uses Armadillo. I get very close but not identical results: Mathematica MachinePrecision: {1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 1.98037*10^13, ...

0

From FindMinum's point of view the objective function you are minimizing is a "black-box function", meaning it is an oracle, that yields values, but not derivatives, or Hessians. If the method chosen by FindMinimum to solve the problem requires a gradient of the objective function at a point, finite differences are used to approximate it. The multitude ...

4

ClearAll["Global`*"] T = Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, 3}]; Part[Part[Part[T, 1], 1], 1] = E^(3 (KK + h)); Part[Part[Part[T, 1], 3], 3] = 3 E^(-KK + h); Part[Part[Part[T, 2], 2], 2] = E^(3 (KK - h)); Part[Part[Part[T, 2], 3], 3] = 3 E^-(KK + h); T = Normal[Symmetrize[T]]; h = 0; KK = 10; im = Partition[Flatten[TensorContract[TensorProduct[T, T], ...

2

As I noted in a comment above, it is quite possible for the solution of a nonlinear PDE to become singular at finite t, and that appears to be occurring here. And, as noted by Dr. Wolfgang Hintze, the right side of the PDE, when integrated over {x, -6, 6} is zero. So, the integral of f must be a constant, and indeed it is. For the parameters given in the ...

10

EDIT #2 My error was useful. It brought me to the conclusion that the difficulties in solving the PDE of the OP are due to the drift term $$\frac{\partial (x u(x,t))}{\partial x}$$ If the drift term is included, many boundary problems are ill defined. It turns out that there are cases where mathematically there is only a trivial solution u = 0 but ...

13

Numerics in Mathematica can be as precise as you like. However, precision comes at price; you pay for it in computation time and in additional coding effort. In Mathematica there are several computational classes of non-complex numbers, which form a tree like this. The computation you made was made with machine reals because you included 0.5 as a term. ...

0

Today I believe I encountered the same problem when trying to reproduce the result of this paper about Lamb's problem and solutions mentioned above doesn't help in my specific case. After struggling for a while I managed to find a work-around and I think it's worth sharing. In short, if the integrate contains singular point(s) in addition, you may need to ...

1

I think that the function f in Mr. Wizard's answer linked above may be what you want. But let me ask you, when you say you want N[1/Sqrt[3], 1] which is 0.6, do you want that precision maintained throughout calculations? That is, if I enter 1/.6 it returns 1.66667 with a repeating decimal representation. But if I force the precision to be 1 on this ...

Top 50 recent answers are included