Hot answers tagged numerical-integration
22
There is an (undocumented?) feature of NDSolve which is handy for exactly this purpose: You can add more than just the start and end of the integration interval and enforce that these points will be met. The result is like you would run NDSolve on each of the corresponding intervals with the starting conditions given by the end point of the previous ...
19
Perhaps setting the difference order to
"DifferenceOrder" -> "Pseudospectral" is what you are looking for:
showStatus[status_] :=
LinkWrite[$ParentLink,
SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],
ToString[status]]];
clearStatus[] := showStatus[""];
clearStatus[]
nxy = 33;
sol = NDSolve[{D[\[Rho]g[t, x], t] + D[\[Rho]g[t, x] u[t, x], ...
17
Generally speaking, you can recognize a list because it'll have List as its Head. For example:
Head[{1,2,3}]
will return List.
For your example conditional where you want to change what you do based on the Head of the resulting expression, you can use Switch, such as in:
Switch[result,
_List, what you want to do with a list,
_, what you ...
17
My variant of Szabolcs code. It doesn't need an extra package:
sol = First[
NDSolve[eqns, {a, b}, {t, 0, 1000}, Method -> {"EventLocator",
"Event" -> Abs[a'[t]] +Abs[b'[t]] < 10^-5,
"EventAction" :> Throw[end = t, "StopIntegration"]}]];
Plot[Evaluate[{a[t], b[t]} /. sol], {t, 0, end}]
As you can see it makes use of the ...
16
You can use the EventLocator method of NDSolve.
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
eqns = {Derivative[1][a][t] == -a[t] - 0.2` a[t]^2 + 2.1` b[t],
Derivative[1][b][t] == a[t] + 0.1` a[t]^2 - 1.1` b[t],
a[0] == 0.5`, b[0] == 0.5`};
sol = First@
NDSolve[eqns, {a, b}, {t, 0, 1000},
Method -> {"EventLocator",
...
16
You are trying to solve an elliptic, purely spatial problem. As described in the documentation on numerical PDEs, Mathematica uses the numerical method of lines, which requires a temporal variable with initial (not boundary) conditions. That is the proximate cause of your problem.
Honestly, in this case, you might consider trying another system that ...
16
Time-dependent case
in the time-dependent case, $[H(t),H(t')]\neq0$ in general and we need to time-order, ie, the operator taking a state from $t=0$ to $t=\tau$ is $U(0,\tau)=\mathcal{T}\exp(-i\int_0^\tau dt\, H(t))$ with $\mathcal{T}$ the time-ordering operator. In practice we just split the time interval into lots of small pieces (basically using the ...
15
You can always separate your inner integrals, convert them to functions and use in NIntegrate:
i1[z_?NumericQ] := i1[z] = NIntegrate[-y, {y, 0, z}]
i2[x_?NumericQ] := i2[x] = NIntegrate[Exp[i1[z]], {z, -∞, x}]
NIntegrate[x i2[x], {x, -5., 5}]
(* 30.0795 *)
14
NIntegrate performs a certain symbolic processing of the integrand to detect discontinuities, singularities, to determine the method to choose and so on.
If you know the integrand pretty well, the way to reduce the overhead is to set the method explicitly, set its SymbolicProcessing suboption to 0 (to allow to time spent on the preprocessing), and to add ...
14
As it turns out, the designers of NDSolve[] have precisely anticipated this sort of use; this is where you can use the NDSolve`StateData framework.
To use acl's example:
(* prepare PDE *)
state = First[NDSolve`ProcessEquations[{D[u[t, x], t] == D[u[t, x], x, x],
u[0, x] == 0, u[t, 0] == Sin[t], u[t, 5] == 0}, u, t, {x, 0, 5}]];
(* go up to t = 2 *)
...
13
There is the function NFourierTransform[] (as well as NInverseFourierTransform[]) implemented in the package FourierSeries`. The function, as with the related kernel functions, takes a FourierParameters option so you can adjust computations to your preferred normalization as needed. For your specific normalization, you apparently want the setting ...
13
Some frames from my version of the animation:
Here's the code I used:
orbit[posStart_?VectorQ, derStart_?VectorQ] :=
Block[{c = -Rationalize[6.672*^-11*7*^17], x, y, z, t},
{x, y, z} /. First @ NDSolve[
Join[Thread[{x''[t], y''[t], z''[t]} ==
c {x[t], y[t], z[t]}/Norm[{x[t], y[t], z[t]}]^3],
...
12
This approach finds equilibrium by checking that all derivatives up to the order of the differential equation are below a threshold. Following the template (defined below) suggested by the OP, here is an example for a damped harmonic oscillator:
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
eqns1 = {a''[t] == Pi^2/2500 - (Pi^2*a[t])/2500 - ...
12
If you know the equation defining your ellipsoid you could use Boole[] to constrain the integration domain :
myF[x_,y_]=Abs[x+y]
NIntegrate[Boole[(x/3)^2 + (y/2)^2 <= 1] myF[x,y], {x, -5, 5}, {y, -5, 5}]
Note that this will actually prevent myF[x, y] from being evaluated outside the domain specified by Boole. This feature of NIntegrate is described ...
12
noeckel’s answer on StackOverflow is spot on. This is not a Mathematica issue, this is a mathematical issue. Namely, Mathematica is giving you the correct solution to the system of differential equation and boundary conditions given. The conditions given (and in particular the derivative imposed at the origin) are incompatible with the expected decay. Bear ...
12
This is fixed in version 9.
This came up on MathGroup before. Since it hasn't been fixed for so long, I wasn't sure if it was really a bug, so I did some spelunking (and some speculation) today to find out what's happening. To jump to the end: I think it's a bug.
First, let's see what arguments does LogLinearPlot really pass to the function:
...
12
You can modify the global system variable $Assumptions, to get the effect you want:
$Assumptions = aa[t] > 0
Then
Integrate[D[yy[x, t], t]^2, {x, 0, 18}]
10.1601 Derivative[1][aa][t]^2
This may, however, be somewhat error-prone. Here is how I'd do this with local environments. This is a generator for a local environment:
...
12
As the previous answer shows while dealing with this type of problem in Mathematica one must use _?NumericQ in the function definition. Once you define the function as Mr.Wizard♦ has instructed, it is pretty straightforward to call NMinimize or FindMinimum. I would also like to scale your function fcc by $10^6$. However FindMinimum is better suited for this ...
12
Here is a code that is about 2 orders of magnitude faster. We will use a finite element method to solve the issue at hand. Before we start, note however, that the transition between the Dirichlet values should be smooth.
We use the finite element method because that works for general domains and some meshing utilities exist here and in the links there in. ...
11
David's method is the route I would go if I didn't want to use Quantile for some reason.
Quantile[StudentTDistribution[49], .95]
==> 1.67655
If you really want to use PDF's to demonstrate this you can still use FindRoot by first setting the integration up so that it only accepts numeric input.
f[y_?NumericQ] := ...
11
I think you intended to use {li, 200, 800} instead of {li, 800, 200}.
If you do so, then you could visualize the result :
ListLinePlot@dnFpoints
Moreover I would rather define daF in the following form :
daF[l_]:= 500 * 0.28 Exp[-((l - 500)/90)^2]
c = 3 10^8;
Edit
Instead of using Table of dnFpoints I add an alternative method for calculation of ...
11
Why not use elliptical coordinates? That's what they are there for.
For example, if your function is f[x_,y_], then you define the coordinate transformation
x[u_, v_]:= Cos[v] Cosh[u];
y[u_, v_]:= Sin[v] Sinh[u];
The Jacobian is Sin[v]^2 + Sinh[u]^2, so you simply do the integral like this:
NIntegrate[
f[x[u, v], y[u, v]] (Sin[v]^2 + Sinh[u]^2), {u, ...
11
For total control over integration, sweep over the ellipse along the eigendirections. To get the integral right, it is essential to use unit-length eigenvectors. This code shows a discrete version of the sweep, to help you visualize it:
m = {{2, 1}, {1, 3}}; (* Matrix *)
c = {3, 5}; (* Center *)
ev = Sqrt[Eigenvalues[m]]; (* Semi-axes *)
evec = ...
11
You can use FindRoot :
func[T_?NumericQ] := NIntegrate[1/(E^(1/(\[Lambda] T)) - 1), {\[Lambda], 200, 220}]
sol = FindRoot[func[T] == 1000, {T, 0.1}]
(* {T -> 0.240468} *)
func[sol[[1, 2]]]
(* 1000. *)
FindRoot needs an initial guess for the unknown; you can get a rough estimate by plotting func for instance.
11
This is how you manually invoke "LevinRule" when you know part of the integrand is a rapidly oscillatory function satisfying a linear ODE:
First, a rapidly oscillatory function:
In[25]:= osc =
y /. NDSolve[{y''[x] - (x^2 - 3 x) y'[x] + 10000 y[x] == 0,
y[0] == 3, y'[0] == 1}, y, {x, 0, 5}] // First
Out[25]= ...
11
You can define function to only evaluate on explicit reals:
F[x_Real] := NIntegrate[f[y], {y, 0, x}];
Plot[F[x], {x, 0, 100}]
But a more efficient way would be to use NDSolve instead of NIntegrate:
F[x_] = F1[x] /. NDSolve[F1'[x] == f[x] && F1[0] == 0, F1, {x, 0, 100}];
Plot[F[x], {x, 0, 100}]
10
This question is somewhat subjective, but here's my take on it:
The reason the precise methods are mentioned in papers is to make results reproducible. One has to draw a line when it comes to describing methods. Will you mention what method you used to add or multiply numbers on a computer? What if the numbers are huge and you used FFT-accelerated ...
10
What I'd do is to define the solution as a regular function :
sol[t_] = NDSolve[{Derivative[2][y][t] + 0.1*Derivative[1][y][t] + Sin[y[t]] ==
0, Derivative[1][y][0] == 0, y[0] == 1}, y[t], {t, 0, 10}, Method -> {"EventLocator",
"Event" -> Derivative[1][y][t] - y[t]}][[1, 1, 2]]
Now you can use it as any other function :
sol[2.4985352432136567]
...
10
The issue is that event detection itself adds plenty of overhead. We can see this by comparing this:
model = {x'[t] == x[t] (1 - x[t]) - x[t] y[t],
y'[t] == x[t] y[t] - y[t], x[0] == 0.5, y[0] == 0.5};
perturb = WhenEvent[Mod[t, 1], {x[t] -> x[t], y[t] -> y[t]}];
eventPoints = Reap[NDSolve[{model, perturb}, {x, y}, {t, 0, 10000},
StepMonitor ...
10
You can get the curve in polynomial implicit form as below.
poly =
GroebnerBasis[{x^2 - ct, y^2 - st, ct^2 + st^2 - 1}, {x, y}, {ct,
st}][[1]]
(* Out[290]= -1 + x^4 + y^4 *)
To get the area, integrate the characteristic function for the interior of the region. That that's where the polynomial is nonpositive (just notice that it is negative at the ...
Only top voted, non community-wiki answers of a minimum length are eligible
