New answers tagged

1

As mentioned in the comments that DSolve will be unable to solve this nonlinear coupled system of two ODE's. Instead you can use NDSolve. Here is my try omega = 1.83465945; a0 = 0; epsilon = 0.5; Eq1 = y1''[x] == (Sqrt[(1 - epsilon^2)/((Sqrt[1 + a0^2] + y2[x])^2 - epsilon^2*(1 + y1[x]^2))]/epsilon^2 - omega^2) y1[x]; Eq2 = y2''[x] == (Sqrt[(1 - ...


1

a bit of an extended comment, but i case anyone doesnt see the issue clearly, LegendreP[46, x] is a 46 order polynomial, with all even powers of x and alternating signs on the coefficients. We can separate out the positive and negative terms: {neg, pos} = { Total@MapIndexed[# x^(4 First@#2 - 4) &, #[[1 ;; ;; 2]]], Total@MapIndexed[# x^(4 First@#2 - ...


11

There appears to be a bug, not in Integrate or in BesselK, but in the vertical-axis Ticks of LogLogPlot. Consider the simple case, LogLogPlot[Exp[x], {x, 10^-10, 1}, PlotRange -> All] as it should be. However, LogLogPlot[Exp[x], {x, 10^-10, 1}, WorkingPrecision -> 50, PlotRange -> All] In fact, any value of WorkingPrecision except ...


3

From the LegendreP help page:


3

If you want to tweak the number of digits displayed in your notebook, run this: SetOptions[EvaluationNotebook[], PrintPrecision -> 10] As noted by Szabolcs, the default setting of PrintPrecision is 6, which is why you're only seeing that many digits in the output, even tho all the digits are still there.


11

Hint: try In[1]:= FullForm[N[1.000001, 10]] which returns Out[1]//FullForm= 1.000001` That tells you the 'rounding' is happening on the front end only, but that the full precision you asked for is still there. Roughly speaking, objects have an internal representation, and the front end 'interprets' this representation to produce the display in a ...


15

This has nothing to do with N. You are observing the fact that by default Mathematica truncates machine precision numbers to 6 digits for displaying them. Enter 1.000001 without N and evaluate it: you'll see the same output (i.e. "1."). You can adjust this in Preferences, Appearance, Numbers, Formatting. The numbers are still stored to full precision, ...


3

I have found my mistake. It was on the initial condition if we take $$\phi(zi)=0$$ and $$\phi'(zi)=0$$ the both methods will give the exact same solution.


1

I Have finally found what was causing the problem. it was with notation. I define my constant a_o but in the equation I have used a_0. Once I changed that it was working.


8

The culprit, as suspected by xslittlegrass, is indeed numerical instability; in particular, this is because of the perverse combination of modified Bessel functions exhibited in the result returned by Mathematica. Using a recurrence identity satisfied by the modified Bessel function of the first kind, we can simplify the expression returned, like so: ...


9

It seems that the analytic result is correct, but the precision is lost when converting it to a number. For example, if we use a higher precision, we get consistent results between numerical and analytical integration: f[a_, b_] = Integrate[x^2 Exp[-a x^2 - b x^4], {x, -∞, ∞}, Assumptions -> {a > 0, b > 0}] g[a_, b_] := NIntegrate[x^2 Exp[-a ...


6

The oddity in this case comes from NSum which is being called in a certain way from NIntegrate. This is a simple example that has roughly the same behavior (note in this case the exact result is known to be $\mp \ln 2$): NSum[(-1)^n/n, {n, 1, Infinity}, Method -> {"AlternatingSigns", Method -> "WynnEpsilon"}, WorkingPrecision -> 32] (* ...


3

Check out the following FullSimplify[(x^2 + y^2) Cos[ 4 ArcTan[y/x]] - ((8 x^4)/(x^2 + y^2) - 7 x^2 + y^2)] (* 0 *) The two expressions are identical so it is not suprising that they produce identical plots.


18

There are several important things about the way computer systems represent real numbers, which most of the time can be blithely ignored, just like the safety of bridges in the United States. One important thing is that numbers are discrete. With regular machine precision (double precision), the mantissa has 53 bits, which provides a lot of resolution. ...


4

Mathematica is checking to make sure that you have enough equations and unknowns, but the way that you've written them out, it only thinks that you have three equations. This is because you've set them up as {5-component vector} == 0, rather than {5-component vector} == {5-component vector}. There are a couple of ways to fix this. One is to use Thread: ...


2

Since your number is quite large, you can use Stirling's approximation to do this. It's also very common to use this approximation in statistical mechanics: For large number $n$ $$\log(n!)\approx n\log(n)-n$$ So in your case $$n=\frac{Nn}{3}$$ then $$ \log\left[\frac{(3n)!}{n!\times n!\times n!}\right]\\ =\log[(3n)!]-3\log[n!]\\ \approx 3n \log(3n)-3n ...


0

It looks like there isn't a good way to bypass Overflow[]. However, for very large factorial calculations, it's useful and incredibly accurate to use Stirling or Nemes approximations, depending on the size of the factorial. @Mathematica devs, an idea- maybe catch overflow errors, tell Factorial to substitute the Stirling or Nemes approximation, then try to ...


7

Spectral methods I present two general ways to approach a second-order linear BVP of the form $$\gamma(x)\, y''(x) + \beta(x)\, y'(x) + \alpha(x)\, y(x) + \varphi(x) = 0,\ y(0) = y_1,\ y(\infty) = y_2$$ By two changes of variables, it can be put into the following forms, including one with homogeneous boundary conditions: $$C(t)\, u''(t) + B(t)\, u'(t) + ...


3

The numerical definitions you give later modify the calculation of the symbolic definitions you use earlier. If you clear all your variables at the start of your calculation, the problem goes away and repeatable results are obtained. In particular, it is the numerical definition of k that seems to muddy the waters. In order to do what you want, a far better ...


3

In this BVP, Mathematica uses the shooting method, and you're shooting at an unstable solution. You can iteratively approach the solution, extending the interval of integration at each step, but at machine precision, you can extend it only so far. I'll set up bc to depend on symbolic min and max, so we can fiddle with them at each iteration: (*Boundary ...


0

The problem that you presented can be solved using FindFit. You are trying to fit the actual function f[x_] := Exp[x]^(1/2) + 2 - Exp[x] + x^5 in the range zero to one with a quadratic. First step is to generate the data. Creating it within the range 0 to 1 satisfies the fist constraint. data = Table[{x, f[x]}, {x, 0, 1, 0.02}]; ListPlot[data] ...


1

Oh! I'm coming to late! r1 = 1.8; r2 = 2.2; ode1 = h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r]); ode2 = p'[ r] == -(2 p[r]/r) - (h'[r]/h[r])* p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r])); sol = ParametricNDSolve[{ode1, ode2, m[r] == (r/2)*(1 - h[r]), m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}] One way to find ...


2

r1 = 1.8; r2 = 2.2; ode1 = {h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r])}; ode2 = {p'[r] == -(2 p[r]/r) - (h'[r]/h[r])* p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]))}; hpm = ParametricNDSolveValue[{ode1, ode2, m[r] == (r/2)*(1 - h[r]), m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}] Plot[Evaluate[Through[hpm[0.8, 1.][t]]], ...


10

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


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


5

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


7

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 /@ % ...



Top 50 recent answers are included