New answers tagged

7

I thought this was described in some NDSolve tutorial, but I can't find it. It is mentioned in Only final result from NDSolve, which @andre314 linked in a comment. The call NDSolve`ProcessEquations[{y’[x] == y[x], y[a] == 1}, y, {x, b, c}] with an initial condition at x == a and interval b <= x <= c will save only the solution between x == b and x == ...


4

We should note first that point $X=0$ is a branching point for the first order equation. Therefore we need to make special step before RK4. Using initial condition we calculate, that $C=D=0$, therefore first order equation looks like $$y'=\pm \sqrt {2 y^3+y^2}$$ Code for RK4 we can organize with Module as follows rk4[f_, h_, x0_, y0_, n_] := Module[{k1, k2,...


5

Changing the numerical integration strategy from the default "GlobalAdaptive" to "LocalAdaptive" gives the expected behavior (and results much closer to the expected ones): NIntegrate[Exp[-(x - y)^2], {x, -100, 100}, {y, -100, 100}, Method -> "LocalAdaptive"] (* 352.491 *) NIntegrate[Exp[-(x - y)^2], {x, -200, 200}, {y, -200,...


0

This integral can be don analytically. Then you do not have to worry about precision. And you see that the exponential term can be safely neglected. Integrate[Exp[-(x - y)^2], {x, -200, 200}, {y, -200, 200}]


0

First look at the function you want to integrate Plot3D[Exp[-(x - y)^2], {x, -5, 5}, {y, -5, 5}, AxesLabel -> {x, y}] The plot shows integral must depend on the range of integration! That's why you cannot expect unique solution.


5

A similar but lower level approach to what Tim suggested is to use DiffisionPDETerm and ConvectionPDETerm vars = {u[t, x, y], t, {x, y}}; eqn2 = D[u[t, x, y], t] + DiffusionPDETerm[vars, 0.1] + ConvectionPDETerm[vars, {y, -x}] == 0; Boundary conditions are then: bc = DirichletCondition[u[t, x, y] == Sin[Pi*x*y], True] The rest is then the same: (* ic =...


6

As I alluded to in the comments, I described a complete description of setting up a well-formed PDE system for a transient convection-diffusion problem in my answer here 237643. The workflow for rectangular domain using the HeatTransferPDEComponent is shown in the following: vars = {u[t, x, y], t, {x, y}}; pars = <|"ThermalConductivity" -> 0....


0

Update Your piecewise definition for YO is improperly written. However, I would recommend not using Piecewise at all. I would rewrite your code like this: θ = 0; λ = 1.2398/(19.5 10^3); f = 4.72 10^3; Subscript[δ, 1] = 1.274/10^6; Subscript[δ, 2] = 4.304/10^6; Subscript[bt, 1] = 5.254/10^9; Subscript[bt, 2] = 2.435/10^7; χ1 = -2 Subscript[δ, 1] + 2 I ...


2

About the largest value of inf that yields a result for the code as given in the question is 0.6. SetOptions[Plot, ImageSize -> Large, LabelStyle -> {15, Bold, Black}]; inf = 3/5; wp = 30; sol = NDSolveValue[{test1, test2, x[inf] == 2, y[inf] == 4, x[-inf] == 1, y[-inf] == 2}, {x[t], y[t]}, {t, -inf, inf}, WorkingPrecision -> wp]; g = D[...


2

$Version (* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *) Clear["Global`*"] Use Piecewise rather than If for numeric functions Function1[x_, y_] = Piecewise[{{Exp[-x^2 + y^2], 4 < x^3 + y^4 - x^2 < 5}}]; int1 = NIntegrate[Function1[x, y], {x, 0, 100}, {y, 0, 100}] (* 0.597656 *) reg = ImplicitRegion[ 4 < x^3 +...


3

With the data: Function1[x_, y_] = If[4 < x^3 + y^4 - x^2 < 5, Exp[-x^2 + y^2], 0]; DataInt = Flatten[ParallelTable[{x, y, Function1[x, y]}, {x, 0, 100, 0.05}, {y, 0, 100, 0.05}], {2, 1}]; FunctionInt[x_, y_] = Interpolation[DataInt, InterpolationOrder -> 1][x, y]; You could try: ir = ImplicitRegion[{4 < x^3 + y^4 - x^2 < 5, 0 < x ...


4

You can use Integrate to accumulate an InterpolatingFunction: Integrate[FunctionInt[x, y], x, y] /. {x -> 100, y -> 100} (* 0.601376 *) Addendum: If the interpolation grid is regular as in the example (spacing = 0.05), here's a quick way using a manual trapezoidal rule: With[{fvals = Head[FunctionInt[x, y]]["ValuesOnGrid"]}, Nest[Total[...


3

Find the domain over which the values are not zero: minmaxes = MinMax /@ Transpose@Select[DataInt, #[[3]] != 0 &] (* Out: {{0., 2.1}, {0., 1.5}, {0.0121552, 8.67114}} *) Select a square region that includes that domain and interpolate over that: squareint = Interpolation[ Select[DataInt, 0 <= #[[1]] <= 2.1 && 0 <= #[[2]] <= ...


3

The main part of your integrand is Sin[t]^40000. This produces an infinity of very sharp peaks. Here is part of the graph: I doubt if there is a purely numerical algorithm, that can integrate this from 0 to infinity. Does this integrand appear in a real problem or is it rather a made up contraption to fool numerical routines. In the first case, some ...


2

thisstep = 0; laststep = 0; stepsize = 0; First@NDSolve[{eqns, Y1[0] == X0, Y1'[0] == X0d}, Y1[t], {t, 0, tmax}, MaxStepFraction -> 1/15, StepMonitor :> (laststep = thisstep; thisstep = t; stepsize = thisstep - laststep;), Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}, ...


1

12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020) gives these results using antiderivate with integration limits 0 and Pi: Plot[{Re[NIntegrate[Sqrt[(37 - 45*37*(x^2/(74*150)))^2*Sin[t]^2 - (40 - 27*37*(x^2/(16*150)))^2*Cos[t]^2], {t, 0, Pi}, WorkingPrecision -> 50]], Im[(1/400)*(32000 - 333*x^2)* EllipticE[(1900160000 - 28416000*x^2 + ...


8

If the DQM is implemented correctly, then this may be an essential limitation of the method. I knew nothing about DQM, but scanning this paper, I have the feeling the method is similar to Pseudospectral. Indeed, a quick test shows that, the weighting coefficient matrix of 1st order derivative in DQM is exactly the same as the differentiation matrix of 1st ...


10

This is an interesting example where new improvements in version Mathematica 12.2 of elliptic functions handling appear important to get a correct result. TraditionalForm[ intd[x_] = FullSimplify[ Integrate[ Sqrt[(37 - (45 37 x^2)/(74 150))^2 Sin[t]^2 - (40 - (27 37 x^2)/(16 150))^2 Cos[t]^2], ...


8

Asymptotics (integration by parts) Since the integrand is highly oscillatory for large $s$, we can obtain an asymptotic expression by performing integration by parts. Observe that $$ \begin{align} & \int_0^1 x^2 \mathrm{e}^{\mathrm{i} f(x) s} \,\mathrm{d}x \\ &=\int_0^1 \frac{x^2}{\mathrm{i} f'(x) s} \,\mathrm{d}\left( \mathrm{e}^{\mathrm{i} f(x) s} \...


4

Once again demonstrating the power of Hadamard's maxim, let us use a parabolic contour to evaluate this oscillatory integral: parabolic[a_, x_] = Simplify[InterpolatingPolynomial[{{0, 0}, {1/2, -a}, {1, 0}}, x]]; With[{s = 1*^7, a = 1/2}, NIntegrate[With[{x = x + I parabolic[a, x]}, x^2 Exp[I Log[1 + x]/x s]] (1 + ...


1

I assume that x[t] in the definition of L2 is a typo and should read x[r]. With this assumption, you can simply define L2 as: L2 = r^3 ((V2[x[r], r])*(1 + ((1/2)*(x'[r]^2))) + (1/Pi)*(Pi*(800) - x[r] + (1/2) Sin[2*x[r]])) The equations then read: eq={D[D[L2, x'[r]], r] == D[L2, x[r]], x[0] == 0, x'[0] == 0} Furthermore, the syntax of NDSolveis ...


6

Here are two ways, a simple change of variables and an assist to the Levin Rule. Change of variables Block[{f, s}, integrand = ComplexExpand[ x^2 Exp[I f[x] s], TargetFunctions -> {Re, Im} ] Dt[x]; f[x_] := Log[1 + x]/x; s = 10^7; NIntegrate[integrand /. x -> 1/u /. Dt[u] -> 1, {u, Infinity, 1}] ] (* -2.48219*10^-9 + 5.17734*10^-7 I *...


0

If you use the method "LocalAdaptive", the integration routine of MMA will put more sample points in the region where the integrand oscillates strongly. Here is an example: NIntegrate[x^2 Exp[I Log[1 + x]/x 10000000], {x, 0, 1}, Method -> "LocalAdaptive"] This gives: 9.90376*10^-8 + 0.0000697064 I


6

We can draw y=f[x] by ParametricPlot[{x, f[x]}, {x, 0, 2}] and draw its inverse x=f[y] by ParametricPlot[{f[y], y}, {y, 0, 2}] f[x_] := 2 A ArcTanh[(# A)/Sqrt[-1 + #^2 B]] + A Log[1 + #^2 (A^2 - B)] - 2 Log[# B + Sqrt[-1 + #^2 B]] &[ x] /. {A -> 0.2, B -> 0.3} // Abs; ParametricPlot[{{x, f[x]}, {f[x], x}}, {x, 0, 2}] We can also ...


12

The problem is a combination of classic stiffness and weak singularities due to linear interpolation of the data. The stiffness needs a stiff solver plus a higher working precision, as @Ulrich has observed previously and we will explain below. For the weak singularites, NIntegrate has a method "InterpolationPointsSubdivision" to deal with ...


6

What you call g[T] in your introduction seems to be equal to g[t]==-rate[t]/(1 + 2.2*10^-7 t^6/mass^2)^2 jacobian[t]. This values are of order 10^26! You need to increase accuracy inside NDSolve dramatically. Try F = NDSolveValue[{D[f[T], T] jacobian[T]^-1 == fDynamics[f[T], T], f[Tini] == 0}, f, {T, Tfinal, Tini}, WorkingPrecision -> 50] //Quiet ...


3

Use Cartesion coordinate and define a region to cut the surface and cut the domain of integrate. f[x_, y_] := BesselJ[0, r]^2 /. r -> Sqrt[x^2 + y^2]; plot = Plot3D[f[x, y], {x, y} ∈ Disk[{0, 0}, 12], PlotPoints -> 50, PlotRange -> All]; reg = ImplicitRegion[x <= -5 && Sqrt[x^2 + y^2] <= 12, {x, y}]; reg // Region NIntegrate[f[x, ...


7

Since this code is implementation DQM for cantilever beam then we need to put right boundary condition to make this code stable with number of grid points Np changing. This is small modification only but it is works for any Np, for example Np = 20; G1 = 0.05; Ω = 1; μ = 0.05; tmax = 10; a = 30; ii = Range[1, Np]; X = 0.5 (1 - Cos[(ii - 1)/(Np - 1) π]); Xi[...


1

Small round-off errors in t2[l, w] and t1[l, w] lead them to have small imaginary parts, randomly changing as w changes. NIntegrate will integrate between them along a line in the complex plane. In turn, the small, random imaginary parts are propagated through the evaluation of q[r, l, w] = Sqrt[...]. For certain values of r the argument of the square-root ...


1

A generation ago Plot used to let all (or many) messages escape, so that if you plotted a function like 1/x, you wouldn't be surprised at a divide-by-zero error message. You got used to typing Plot[1/x, {x, -1, 1}] // Quiet. Nowadays, Plot is very aggressive at suppressing message, too aggressive in this case, imo. Modifying @Szabolcs's code, here's a way to ...


4

I will use the definition of g given in the answer to question you link to. g[t_, x_] := t^3 - t + x^2 f[t_?NumericQ] := NIntegrate[g[t, x], {x, -1, 1}] minPt = With[{min = FindMinimum[f[t], t ∈ Interval[{0, 1}]]}, {min[[2, 1, 2]], min[[1]]}] {0.57735, -0.103134} Plot[f[t], {t, 0, 1}, Epilog -> {Red, AbsolutePointSize[8], Point @ minPt}] This ...


6

I think that ODE systems with $\| X'\| \sim O(\| X \|^2)$ tend to be unstable, that is, a small rounding error has a chance to cause a solution to blow up. First, boundary-value problems (BVPs) are not guaranteed to have solutions, and without a proof or evidence that a solution exists, difficulty in solving one should raise the question whether there is a ...


5

We can change the variable 2h-h^2->z by hand. Reduce[z == 2 h - h^2 && 1 > h > 0, h, Reals] 0 < z < 1 && h == 1 - Sqrt[1 - z] 1/D[2 h - h^2, h] /. h -> 1 - Sqrt[1 - z] // Simplify 1/(2 Sqrt[1 - z]) It means that the change variable is one to one and $0<z<1$, $\mathrm{d}h=\frac{1}{2 \sqrt{1-z}}$ Integrate[ ...


10

Let NDSolve do the job one time (instead of Integrate many times) Edit Make it faster. AbsouteTiming for all calculations == 0.015 seconds. (tpw[t_] = Total@pulses; gsol = g /. First@NDSolve[{g'[t] == tpw[t], g[0] == 0}, g, {t, 0, 250}]; norm = gsol[250] ;(*11.*) {t05 = t /. First@FindRoot[gsol[t] == .05 norm, {t, 0, 150}], t95 = t /. First@...


8

Should be a bit faster... Long story short: Supply the Jacobian of your equation to FindRoot, too. ϕ = t \[Function] Evaluate[Total@pulses]; ClearAll[Φ]; Φ[x_?NumericQ] := NIntegrate[ϕ[t], {t, 0., x}, PrecisionGoal -> 8]; rhs = 0.05 Φ[tmax]; (*An approximate inverse of Φ*) tlist = Subdivide[0., tmax, 100]; Ψ = Quiet[ Interpolation[ Transpose[...


2

Integrate may use transformations that are not numerically stable. In this case, it uses Simplify, which I suspect factors out the constant term in the exponential, $e^{k+ax}=e^ke^{ax}$. Actually, Simplify seems to do this in a clever way, $c\,e^{-k+ax}=(c/e^k)\,e^{ax}$. The $e^k$ factor overflows and is automatically converted to arbitrary precision. I ...


4

This is an extended comment versus an answer. Following @user21's link, it suggests refining the mesh at the boundaries. The features are very sharp at the Y = 0 boundary. I am going to apply anisotropic meshing and attempt to capture the features. Furthermore, I'm going to shrink the y-dimension since the solution falls off very quickly. Helper functions To ...


1

Let us first simplify the expression under the integral: expr[j_,x_]:=(a[[j]]/x)*((x/a[[j]]) + (x*c[x][[j]]) + (x^2/2)*c[x][[j]]^2)* BesselK[1, Sqrt[x]] // Simplify; And now let us evaluate the integrals AbsoluteTiming[ NIntegrate[expr[1, x], {x, 0, \[Infinity]}, MinRecursion -> 4, AccuracyGoal -> 5]] (* {0.06447, 200011.} *) AbsoluteTiming[ ...


Top 50 recent answers are included