New answers tagged

0

There's not enough information in the question, but this approach fits one interpretation of the description of the situation. One can do the integral in the original NDSolve calls that generated rad and mm. (This may be what J.M. was suggesting in a comment. Or he may have meant executing a separate NDSolve call. I would recommend the approach given ...


0

If you solve for the second derivatives, you won't have to use "EquationSimplification" -> "Residual" and things will work ok. Solving for the second derivatives be faster if you start with exact coefficients. Also, if you clear l, solve for the derivatives, and then substitute a value for l, Solve won't choke on the algebra. The long time it takes is ...


2

Making my comment an answer: Try to extract the mesh from the eigenfunction and use that for integration. Something like: NIntegrate[ circlepluck[x, y]*circlefuncs[[i]][x, y], {x, y} \[Element] circlefuncs[[1]]["ElementMesh"]] This would switch off the adaptive mesh refinement.


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


2

It looks like the best approach is to use spherical coordinates: ω1 = Sqrt[m1 + r^2]; ω2 = Sqrt[m2 + r^2]; f = 1/(ω1 ω2) 1/(ω1 + ω2 + e) + 1/(2 ω1 ω2) 1/(ω1 - ω2 - e) + 1/(2 ω1 ω2) 1/(ω2 - ω1 - e) (* ==> 1/( 2 Sqrt[m1 + r^2] Sqrt[ m2 + r^2] (-e + Sqrt[m1 + r^2] - Sqrt[m2 + r^2])) + 1/( 2 Sqrt[m1 + r^2] Sqrt[ m2 + r^2] (-e - Sqrt[m1 + r^2] + ...


1

Your problem is twofold. The first part is very subtle and sneaky. Your "c" characters are actually typed in as the Unicode Character 'cyrillic small letter es' (U+0441) in the NDSolve expression, as you can see from the following: Inactive[NDSolve][{ I*с00'[t] == с10[t] + с01[t] + с00[t] + с00[t], I*с01'[t] == с11[t] + с00[t] + с01[t] - с01[t], ...


0

I think this integral can be done analytically. I have ignored the 0.001 term as I suspect this is merely intended to improve the numerical behaviour. Assuming[a > 0, Integrate[w E^(-w/a) Sin[(w) 100000]/((w)^2), {w, 0, Infinity}]] gives the value ArcTan[100000 a]


0

Assuming principal value: Block[{a = 0.001, w0 = 0.001}, NIntegrate[ w E^(-w/a) Sin[(w - w0) 100000]/((w - w0)^2), {w, 0, w0, ∞}, MaxRecursion -> 20, AccuracyGoal -> 10, Method -> "PrincipalValue"] ] (* -0.000052312 *) Check (seems stable with increased precision): Block[{a = 0.001 // Rationalize, w0 = 0.001 // Rationalize}, ...


5

Maybe commandeer FindRoot with the Villegas-Gayley trick: Updated, with the order of the steps taken by FindRoot saved in icsteps. The results of FindRoot, as saved by NDSolve and shown below as DownValues[], have been sorted by Mathematica and are not in the order in which there were called. This update stores the order in icsteps. Clear[x, t]; Internal`...


4

I think this gives what you want: We construct the integral inside the definition of H by adding another ODE to the NDSolve system, which I called logH. This in fact calculates the integral from ic, not from 0. So to define H we need to subtract logH[0] from logH[t] before exponentiating. This should be a much more accurate (and faster) way of computing ...


2

Here is a function for multiplying InterpolatingFunctions generated from a single NDSolve (so that the coordinate grids are the same, as well as one-dimensional). This yields a single InterpolatingFunction that interpolates the product. It carries over derivative information, too. (I've done this before on the site, I think.) Anyway, Integrate on an ...


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


3

You have many superfluous sets of {} that generate unexpected output in your code. In particular, the Interpolation functions generated by NDSolve were not returning a scalar value, but instead a unidimensional vector, i.e. a list containing a single value instead. That was probably an unintended consequence of the extra sets of braces in the definitions of ...


3

Clear["Global`*"] yy = 10^-4; rr = 0.999; xx = 10^-15; zz = 10^-4; mm = 10^-4; yy + rr + xx^2 + zz - mm^2 - zz^2/24 ic = -17.5 s = NDSolve[{D[y[t], t] == (3 y[t])/5 - (12 m[t]^2 y[t])/5 + (2 r[t] y[t])/ 5 - (6 x[t]^2 y[t])/5 + (3 y[t]^2)/5 + (7 y[t] z[t])/ 5 - (y[t] z[t]^2)/10, D[r[t], t] == -((2 r[t])/5) - (12 m[t]^2 r[t])/5 + (2 r[...


0

As mentioned in the comments above, your understanding for function is wrong. You'd better put more effort in learning the core language of Mathematica, or even consider revisiting the definition of function in math. The following is a not-the-simplest but correct way to code your equation set in Mathematica, do read it carefully: Clear["`*"] vei = 1/1000; ...


1

You cannot use square brackets for parentheses and must include initial conditions. eqns = { x'[t] == -z[t] (x[t] + 1.5 Log[z[t]]), y'[t] == z[t] - 1, z'[t] == -(.64) (y[t] + (x[t] + Log[z[t]])/.8), x[0] == y[0] == z[0] == 1} // Rationalize // Simplify; sol = NDSolveValue[eqns, {x[t], y[t], z[t]}, {t, 0, 1000}]; Manipulate[ ...


2

This works Clear[x, y, z] sol = NDSolve[Join[ Thread[{x'[t], y'[t], z'[t]} == {-z[t] x[t] + 1.5 Log[z[t]], z[t] - 1, -(0.64) y[t] + (x[t] + Log[z[t]])/0.8}], Thread[{x[0], y[0], z[0]} == {1, 1, 1}]], {x, y, z}, {t, 0, 100}]; ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 100}]


1

When you see the warning NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help. This may indicate the NDSolve did try to solve this as a purely spatial problem. You can force it so use the MethodOfLines via the option Method -> "MethodOfLines" for time integration. This then solves in ...


1

e = .16; g = .4; w = .97; You must add a initial condition, assume: L[0]=1 sol = NDSolveValue[{y1'[t] == y2[t], y2'[t] == -((2/L[t]) L'[t] + g L[t]) y2[t] - (1/L[t]) Sin[y1[t]], L'[t] == 7 e w Sin[w t + 9 Pi/8]^6* Cos[w t + 9 Pi/8], y1[0] == -1, y2[0] == 1, L[0] == 1}, {y1, y2, L}, {t, 0, 30}]; ParametricPlot[{sol[[1]][t], sol[[2]][t]}, {t, 0, 30}...


2

I have revised the code to use StreamPlot instead of VectorPlot and NDSolve. The latter command created issues when the domain for the solution extended to negative times. I don't know why, but NDSolve produced integral curves that did not follow the direction field produced by VectorPLot as provided by MarcoB's nice reformulation of my original code. I will ...


5

Here is a workaround suggested by the response I received from WRI: {sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., WhenEvent[Norm[{x'[t]}] < 0.0001, {x[t] -> 1.}]}, {x}, {t, 0, 200}, Method -> {"EquationSimplification" -> "Residual"} ] Plot[x[t] /. sol, {t, 0, 200}] Warning: This option works by converting the system to a DAE, ...


3

As I mentioned in comments, your main problem was the placement of the definition of z; that definition does not need to be re-evaluated every time there is a modification within the body of Manipulate since it's already delayed. The correct placement for such a definition would be in the Initialization code for the Manipulate. Similarly, the definitions of ...


1

If we multiply through the second equation by t, NDSolve can find a solution using "EquationSimplification" -> "Residual": sol = NDSolve[{ {x'[t] == y[t], t y'[t] == -y[t] - t Exp[x[t]]}, {x[0] == 1, y[0] == 0}}, {x, y}, {t, 0, 1}, Method -> {"EquationSimplification" -> "Residual"}] Plot[{x[t], y[t]} /. sol // Flatten // Evaluate, {t, 0,...


2

The DE has problems at y[x] == 0 as noted, but also at x == 1, where there is a pole of order 2. This suggests that numerical issues near the singularities could cause substantial error to accumulate. One can use Piecewise to substitute the limiting value at a discontinuity, but this won't cure any numerical instability in the neighborhood of the ...


3

Try: {sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., WhenEvent[Norm[{x'[t]}] < 0.0001, x[t] -> 1.; x[t] -> 1]}, {x}, {t, 0, 200}]


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


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


5

mat[x_?NumberQ] := {{x^2, 1, 0}, {1, x^2, -1}, {0, -1, x}} ei[x_] := Eigensystem[mat[x], 1, Method -> {"Arnoldi", "Criteria" -> "RealPart"}][[1, 1]] NIntegrate[ei[x], {x, 0, 4}] (* 8. *)


3

This answer is just a response to the following part of the question: ……generally, how to speed up the plot involving NIntegrate or is there any principle to do it? When the symbolic processing of NIntegrate isn't necessary for getting the correct result of integration, "SymbolicProcessing" -> 0 is also a general way for speeding up: f[t_] := ...


6

Can NIntegrate remember or make full use of the result of a smaller upper-limit integral? Or generally, how to speed up the plot involving NIntegrate or is there any principle to do it? Below is a solution in the spirit of this request. In short, we find adaptively sampled points, compute integral estimates over the intervals having those points as ...


2

Let us have a single integrand function first: F[x_?NumericQ, xp_?NumericQ] := ( EllipticK[(4 x xp)/(x + xp)^2] Sinh[x] Sinh[ xp])/((x + xp) (Cosh[202] + Cosh[2 x]) (Cosh[202] + Cosh[2 xp])); obtained by expanding f[x, 101.] f[xp, 101.] g[x, xp]. For me (with Mathematica v. 10.3), the integral fails with NIntegrate::inumri: "The integrand (....


0

Replacing NIntegrate with Integrate fixes that specific problem. But then you use symbolic operations again. I don't think this can be done completely with numerical functions.


3

int[a_, b_] = Integrate[Sin[y1[t] y2[t] r]/r, {r, a, b}, Assumptions -> a < b] odes = {y1'[t] == y2[t] y3[t], y2'[t] == -y1[t] y3[t] + int[y2[t], y3[t]], y3'[t] == -0.51 y1[t] y2[t] t, y1[0] == 0, y2[0] == 1, y3[0] == 2} sol = NDSolve[odes, {y1, y2, y3}, {t, 0, 12}] Edit to answer your question in your comment: int[a_?NumericQ, b_?NumericQ,...


5

Actually the observed behavior is in full accord with the HoldAll attribute, just check what happens when there is no such attribute: nIntegrate[x + x, {x, 1, 2}] // Trace {{x + x, 2 x}, nIntegrate[2 x, {x, 1, 2}]} From the above it is seen that the arguments are evaluated before applying the rules associated with the function nIntegrate. The purpose ...


5

From the "Details and Options" section in the docs: NIntegrate first localizes the values of all variables, then evaluates f with the variables being symbolic, and then repeatedly evaluates the result numerically. So I guess this is expected behavior.


7

Because NIntegrate evaluates the integrands before starting the actual integration, in some cases (like this one) it is better to define the integrand function F with the signature F[x_?NumericQ]. BF[n_?NumericQ, x_?NumericQ] := BesselJ[n, x] NIntegrate[BF[9/2, x], {x, 0, 1}] (* 0.000148473 *) Integrate[BesselJ[9/2, x], {x, 0, 1}] %% // N (* Sqrt[2/\[...


2

I tried to make it a comment and then can't control the words. I think the problem is arising from the fact that the value of the function and its derivative is too small to near x=0. Plot[BesselJ[9/2, x], {x, 0, .01}] Plot[Evaluate[D[BesselJ[9/2, x], x]], {x, 0, 0.01}] As you can see the derivative is highly oscillatory at this small range. When ...


3

Expressing the entire business in terms of SphericalBesselJ[] cures the problem: NIntegrate[r With[{x = r BesselJZero[15/2, 1]}, Sqrt[2 x/π] SphericalBesselJ[7, x]]^2, {r, 0, 1/50}] 1.1879560281974252*^-27 The nice thing about SphericalBesselJ[] is that it does not auto-evaluate to a potentially numerically unstable combination of ...


2

EDIT It can get more complicated. Look at a contour given by r[t_] := 1 + 2 Cos[t]; Original post The result of the integral for a (reasonable) closed contour is just 2 I times the area enclosed by the contour. Proof: the integral is $$\int \left(z^*+z\right) \, dz$$ Letting $$z=x+i y$$, $$dz=dx+i dy$$ the integral becomes $$\int 2 x (dx + i dy)=...


2

As requested, comment made into an answer. The problem was with the integration limits in the OP. The points on the circle satisfy (x-x0)^2 + (y-y0)^2 == R^2 so for a point x in the range [x0-R,x0+R] the limits on y are {y0 - Surd[R^2 - (x-x0)^2,2],y0 + Surd[R^2-(x-x0)^2,2]} I use Surd here because MMA may simplify more easily than with Sqrt, since ...


5

First we shall define 'the integration on a curve'. Traditionally, this is defined as integration of f.dl where dl is the length of a small part of the curve. So, using t as a medium, we can explicitly write out the curve's function on a complex plane, here let's assume it's z=2 Exp[I t]. Then we can use t, a real number, as the integration variable, which ...


4

There is some confusion using Plot. Compare the plot of the non-evaluated y in the question with the plot of the evaluated y: Grid[{{Plot[y[0.25, w - 2], {w, 0, 4}, PlotRange -> All, ImageSize -> Medium], Plot[Evaluate@y[0.25, w - 2], {w, 0, 4}, PlotRange -> All, ImageSize -> Medium]}}] We can see that going through the ...


5

A bit speculative, but I think we can see why 171 is a magic number, if we factor the Gamma[n] from the integral: Table[ {n, NIntegrate[Gamma[x] /Gamma[x + n], {x, 1, 2}]} , {n, 168, 173}] {{168, 7.20924*10^-304}, {169, 4.26119*10^-306}, {170, 2.41873*10^-308}, {171, 0.}, {172, 0.}, {173, 0.}} 171 is where the integral fails as you can see the ...


3

Fixing it ... Is there any settings that I can turn on to fix this? Does using larger working precision produce results you expect? F[n_] := NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}, WorkingPrecision -> 60, MinRecursion -> 6, MaxRecursion -> 20] In[10]:= F[170] Out[10]= 0....



Top 50 recent answers are included