# Tag Info

8

Just to contribute to the debate, here is some more evidence that supports the proposition that numerical error is the issue. If we run the integral through various permutations of the ways of making exact and approximate calculations, the pattern I think suggests that numerical error is the reason the OP's integral is so far off. (* the integrand and ...

4

I do not think this is related to floating point errors. I think the closed form solution arrived to in Integrate is not correct. May be wrong branch is taken. To see this more easily, Here is a simpler one (part of the original integral) that gives a symbolic solution, but wrong numerical value for the values when substituted into the expression Clear[a, ...

8

Try not to supply machine numbers to integrals over infinite domains. They can cause errors that build up to the extent you have seen. Either compute the symbolic integral with exact numbers (and then convert it to a numeric value) L = 2; sn = 1; a = 10^(sn/10); b = 10^(sn/10); c = a/100; result = 2*Sqrt[1/Pi]*Integrate[(1/(E^z*Sqrt[z]))*(1 - (a/(a + ...

1

You can use Reduce or CylindricalDecomposition to break down the inequalities into forms that can be used as limits of integration. (Note: One might try BooleanMinimize instead of distr below to distribute And over Or, but it changes the order of the inequalities. Consequently, the intervals of integration end up out of order and NIntegrate does not work. ...

1

One way to proceed is to try out portions and see where the problem lies. For instance: NIntegrate[Boole[0 < e1 < e2 < 1/2 && 0 < e3 < 1/2 && e1 < p < 1 - e1], {p, 0, 1}, {e1, 0, 1/2}, {e2, 0, 1/2}, {e3, 0, 1/2}] 0.0416667 is fine. But the inequalities NIntegrate[Boole[e3 < myV[p, e1, e2] < 1 - e3 < ...

3

Here is another way to solve this issue. I know, it is written everywhere, but this is a common mistake to think that the orifice equation is Cd * Ad * Sqrt[ 2/Rho * ( ps - p1[t] ) ] and could give complex results ! There is no physics going imaginary in the real world. The equation is just wrong. When the pressure reverse, the flow reverse, at least. ...

1

Use a numerical derivative: Clear[band, en, w, fermi, k, T, S, a] << NumericalCalculus k = 86*10^-6; T = 4000; band[en_, w_] := 10000 Exp[-en^2/2 w^2]; fermi[en_, ef_] := 1/(Exp[(en - ef)/(k T)] + 1); S[ef_?NumericQ] := -k NIntegrate[(band[en, w] fermi[en, ef] Log[fermi[en, ef]]) /. w -> 1, {en, -Infinity, Infinity}] ...

1

Ok, so this is my attempt at a short discussion and solution (after getting a little bit wiser). Due to these equations being stiff (Wikipedia) as pointed out by @Nasser, certain numerical methods have difficulty at tracking the solution, and thus, approximations of the orifice equations yield complex solutions at some point when p1[t]>ps and/or ...

4

You could try specifying the event differently (this will trigger when the absolute difference between p1 and p2 is smaller than some value): Ad1=10/1000^2; Ad2=1.5*1000^(-2); Ad3=1.5*1000^(-2); Cd1=0.67; Cd2=0.67; Cd3=0.67; V1=10/1000; V2=10/1000; Rho=875; beta=1000*10^6; ps=100*10^5; Q1=Ad1*Cd1*Sqrt[(2/Rho)*(ps-p1[t])]; ...

4

The problem with WhenEvent has to do with the OP's DE. For an event to be detected, there has to be a point at which the condition is crossed, that is, changes from False for t < t0 to True for t > t0. NDSolve then applies a root-finding algorithm to approximate the value of t0 at which the event occurs. In your DE, the solution p1[t] theoretically ...

0

This is an opportunity for you to clarify your question or give a better example. The coding example that you posted as non-working, does work: tab = {{{1., 3., 7., 5.}, 1.85541}, {{1., 3., 8., 5.}, 1.76612}, {{1., 3., 7., 6.}, 1.99826}, {{1., 3., 8., 6.}, 1.89112}, {{1., 4., 7., 5.}, 0.957483}, {{1., 4., 8., 5.}, 0.868198}, {{1., 4., 7., ...

4

bn[a_, T_, f_] := 2/T Integrate[f Sin[(2 π n)/T t], {t, a, T + a}, Assumptions -> T ∈ Reals] bn[0, L, Piecewise[{{0, 0 < t < L/3}, {w, L/3 < t < 2 L/3}, {0, 2 L/3 < t < L}}]]

1

I've made a number of changes to your code which may speed things up, but I honestly can't say with certainty. Rather than enumerating the changes, I'll just list the code here: ClearAll["Global*"]; integrand[k_?NumericQ, P0_?NumericQ, P1_?NumericQ, rho_?NumericQ, l_?NumericQ] := Module[{x, h}, x = P1*l + P0*(1. - l) - rho; h = HankelH2[0, ...

0

zeta1 = x[t] + I y[t]; eqn1 = (I/zeta1 + 1/Im[zeta1])/(8 Pi Conjugate[zeta1]); s = NDSolve[{x'[t] == Re[eqn1], y'[t] == -Im[eqn1],x[0] == y[0] == 1}, {x, y}, {t, 0, 10}]; Grid[{Plot[# /. s, {t, 0, 10}, Evaluated -> True, PlotLabel -> ##], ParametricPlot[# /. s, {t, 0, 10}, Evaluated -> True, PlotLabel -> ##]} & /@ ...

1

You can try to do it this way: k = Exp[I Pi/4]; f[y_] := y * NIntegrate[ 1/x * HankelH1[1, k*x] * HankelH2[1, x/k], {x, 1, y}] result = NIntegrate[ f[y], {y, 1, 2}, Method -> {Automatic, "SymbolicProcessing" -> None}]

3

2 problems: You were comparing, in the WhenEvent, solution, which had complex value at that t, to real numbers. I used Abs. If this does not work for you, you can use Re, but can't compare complex number to real number using >. Second, your system is stiff, need to use StiffnessSwitching to help NDSolve. d1 = 10/1000^2; Ad1 = 10/1000^2; Ad2 = ...

1

This is the approach I usually use in this situation: f[(a_)?NumericQ, (x_)?NumericQ] := NIntegrate[a*x*k, {k, 0, 1}]; data = {{0, 0}, {0.2, 0.1}, {0.4, 0.2}}; param = FindFit[data, {f[a, x], 2 > a > 1}, a, x] The ?NumericQ predicate will prevent evaluation of f until it gets numerical values, since NIntegrate cannot work with symbolic integrands. ...

1

After correcting some syntax errors and setting consistent boundary conditions: ClearAll["Global`*"]; e = 0.1; pdeset = {Derivative[1, 0][U][t, x] == Derivative[0, 2][U][t, x], Derivative[1, 0][T][t, x] == Derivative[0, 2][T][t, x] + e Derivative[0, 1][U][t, x]^2} ics = {U[0, x] == 0, T[0, x] == 0}; bcs = {U[t, 0] == Sin[t], T[t, 0] == 0, ...

0

The problems were apparently caused by how the parameters were chosen. In the problem at hand they represent physical quantities: parameters starting with r have dimension meter, those with t have dimension second and those starting with v have meter/second. A much more resonable choice would be to use milimeters and miliseconds, i.e. multiplying all ...

5

Observation: First of all I think it is always a useful trick to plot your problematic integrand if possible. It gives us often the clue in case NIntegrate complained about the particular integrand. If we can track down the issue we can often come up with a remedy. Given the following input if one sweeps over the rDet we get the following plots. Given your ...

6

There is a way to get Mathematica to calculate the equation of a transformed cylinder, which can then be used to calculate the volume. First, since you're translating the cylinder, too, I rewrote your transformation to include the translation. We can also define inequalities to define the cylinder. xform[x_, y_, z_, a_, b_] := RotationTransform[a Pi/2, ...

1

Since there is some concern about what the interval of integration should be, why not use a general procedure like this to create all the integrals you need? bIntegrate[f_, {x_, a_, b_}, dx_] := Block[{x0, n}, n = Floor[(b - a)/dx]; x0 = Ceiling[a/dx] dx; Table[Integrate[f, {x, x0 + k dx, x0 + (k + 1) dx}], {k, 0, n - 1}] ] It will split ...

2

Define your scalar product of the h polynomials: p[h[n_], h[m_]] := 16 Pi S Sum[Binomial[8 S + 2, k] (8 S + 1 - n - m - k), {k, 0, 8 S - n - m}] / (2^(8 S + 1) (n + m + 2) (n + m + 1) Binomial[8 S + 2, 8 S - n - m]) Add these properties: p[c_ pol1_h, pol2_] := c p[pol1, pol2] p[pol1_, c_ pol2_h] := c p[pol1, pol2] p[sum_Plus, pol_] := p[#, pol] & /@ ...

1

The following does work in principle but gives a lot of warnings and I don't think you should trust the numbers it gives, although they roughly seem to reproduce the first plot (which is the only I checked but don't know the m it was produced with). I have tried to keep this simple but robust and try to show some techniques that I found helpful in the past. ...

2

This is not the answer to your problem, but perhaps this can help you find a solution to your plotting problem. I would start with defining the equations in terms of the parameters s and m: eq1[m_, s_] = f'''[y] + f[y] f''[y] - m f'[y] - f'[y]^2 - s (f'[y] + y/2 f''[y]) == 0; eq2[m_, s_] = g''[y] - (s/2(3 g[y] + y g'[y]) + 2 g[y] f'[y] - g'[y] f[y]) + ...

0

Your problem is in how you define kx1 and ky1. When you call f1[], x2 and y2 are substituted in the expression they are immediately visible. And since kx1 and ky1 don't explicitly depend on z and x2 and y2, these values aren't substituted. To fix this, you should define kx1 and ky1 as functions: kx1[z_,x2_] = ((2*Pi)/(λ*z))*x2; ky1[z_,y2_] = ...

2

The problem has to do with the sample points used to contruct the plot, and not with NDSolve. Mathematica automatically subdivides segments when the angle is greater than some limit, but it will do so only MaxRecursion times. One can increase PlotPoints to increase the number of initial sample points and increase MaxRecursion to let Mathematica subdivide ...

1

I looked at 5 cases. Conclusion at bottom. ClearAll[x, a, b]; pdf = PDF[NormalDistribution[14, 3.7], x] a0 = pdf[[1]] b0 = pdf[[2, 2, 1]] Case 1 int1 = Integrate[ a Exp[b (-14 + x)^2], x] int1 /. {a -> a0, b -> b0} Simplify[(int1 /. x -> Infinity) - (int1 /. x -> 15)] Case 2 int2 = Integrate[ a Exp[b (-14 + x)^2], ...

1

This is just the working out of the full details of the approach suggested by Michael E2 in his comment of ybelukon's answer. b[i_, j_] := Evaluate @ With[{ f = ((Integrate[(x - i)^i x^j, x] /. Times[coeff_, hyper_Hypergeometric2F1] :> Simplify[coeff, i ∈ Integers && x > i] hyper) /. x ...

1

One can calculate this integral analytically Integrate[(x - i)^i x^j, {x, -5, 5}, Assumptions -> {(i | j) ∈ Integers, j > 0}] ConditionalExpression[(1/(1 + j)) 5^(1 + j) (-i)^i (Hypergeometric2F1[-i, 1 + j, 2 + j, 5/i] + Hypergeometric2F1[-i, 1 + j, 2 + j, -(5/i)] (Cos[j π] + I Sin[j π])), i < -5] Mathematica returns the result for i < ...

Top 50 recent answers are included