Hot answers tagged numerical-integration
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], ...
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 *)
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 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
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
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 ...
9
If you insert a bunch of commands like "Print@First@AbsoluteTiming[...]" in the middle of your function you will see that literally all time is spent on the last line with the Sum[NIntegrate[...]].
To solve your problem insert the following commands into NIntegrate:
Method -> {Automatic, "SymbolicProcessing" -> 0}
Caution: I have personally ...
9
It looks you need to use other NDSolve options to control the spatial grid size to help NDSolve a little.
Trying things, found that by increasing the grid size, now this effect you showed goes away.
Grid size needs to be much smaller when the diffusion coefficient is very small since you are effectively in $\lim{k \to 0}$ now solving $\frac{\partial ...
9
It is possible like so for example:
ClearAll[fn];
fn[z_?NumericQ] := Exp[NIntegrate[-y, {y, 0, z}]];
NIntegrate[x fn[z], {x, -5, 5}, {z, -\[Infinity], x}]
(* 30.0795 *)
but it takes a while to compute.
I used the ability of NIntegrate to handle non-rectangular domains, which is very powerful but not widely known and / or appreciated, it seems. Note ...
9
NDSolve has a slew of options that allow you to control the method. You can find the standard reference here.
There, we learn how to access Euler's method using NDSolve:
Clear[x];
x = x /. First[
NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, x, {t, 0, 10},
StartingStepSize -> 1, Method -> {"FixedStep", Method -> "ExplicitEuler"}]
];
...
9
I hope I see the essence here. You are interested in
the convolution of an interpolated function with a Gauss function
Your underlying data has regular spacings in x-direction and the convolution with a Gaussian is extremely fast implemented in GaussianFilter for discrete data. Why are you making it so complicated when the only thing you have to do is ...
9
I had a play with various Compile options and didn't get anywhere (I managed to make it slower though!). However, you can get a nice little speed boost using ParallelTable. Your original on my machine:
NFourierTransform[f_Function, {kmin_, kmax_}] :=
Interpolation@
Table[{k,
Chop@NIntegrate[f@x E^(-I k x), {x, -Infinity, Infinity}]}, {k,
...
9
You may have some luck playing with various methods. It is also important to decrease the number of integrations - one 4-dimensional integral can often compute much faster than 2 two-dimensional ones (meaning 2d integral of 2d integrals). So, we have:
NIntegrate[Cos[lx]*Cos[qx]+Sin[lx]*Sin[qx],{qx,-Pi,Pi},{qy,-Pi,Pi},
...
9
Let me show how to roll your own numerical solution to a non-linear integral equation using a collocation method. It's fun!
This will involve two approximations. First, we will approximate the function B[x] by its values at n particular points in the range {x, 0, 1}. The integral over x will be replaced by a weighted sum over n, i.e., a quadrature rule. ...
9
You can express your integral in terms of a differential equation and use NDSolve. Since NDSolve builds up the solution as it goes, this is typically much faster.
Clear[y];
y[x_] = y[x] /. First[
NDSolve[{y'[x] == Sin[x], y[0] == 0}, y[x], {x, 0, 10}]
];
t = AbsoluteTime[];
Plot[y[t], {t, 0, 10}]
AbsoluteTime[] - t
8
If you look at the documentation for Precision, it says that if x is the value and dx the "absolute uncertainty", Precision[x] is -Log[10,dx/x]. This, whenever the estimated error is larger than the value, Mathematica will give a negative precision.
Thus, for an estimated error $dx$ and a value $x$ such that $dx/x<1$ here is how the precision as defined ...
8
Here's a slight improvement of b.gatessucks's answer, adding a (mostly effective) initial guess:
te[ii_?NumericQ, opts___] := (\[FormalCapitalT] /. First@FindRoot[
NIntegrate[1/(E^(1/(λ \[FormalCapitalT])) - 1), {λ, 200, 220}] == ii,
{\[FormalCapitalT], (0.5437727672315316 + ii (0.5440978395984463 +
ii ...
8
Check this
f[l_?NumericQ] := Module[
{kk, y, x, out},
kk = NDSolve[{y'[x] == 2*x + 3, y[l] == 0}, y, {x, 0, l}];
out = Evaluate[y[0] /. kk];
out[[1]]
];
res=NMinimize[{f[t], 0 <= t <= 10}, t]
it gives
{-130., {t -> 10.}}
You can check if the answer is correct by plotting your function.
Show[Plot[f[t], {t, 0, 10}, Frame -> True],
...
8
With your definition of a
a[2.4985352432136567]
is Sin[y[t]][2.49854] (which explains the output you get from a[2.4985352432136567] /.sol)
To get Sin[y[2.49854]] you can use
a /. t -> 2.4985352432136567
(* Sin[y[2.49854]] *)
and
a /. t -> 2.4985352432136567 /. sol
gives {-0.556156}.
Alternatively, you can make t an explicit argument ...
8
Had the Euler method not been built-in, one could still use NDSolve[]'s method plug-in framework, which enables NDSolve[] to "know" how to use Euler's method.
Here's how to "teach" NDSolve[] the Euler method:
Euler[]["Step"[rhs_, t_, h_, y_, yp_]] := {h, h yp};
Euler[___]["DifferenceOrder"] := 1;
Euler[___]["StepMode"] := Fixed;
Plugging in the "new" ...
8
As the comments say, you really have to learn correct syntax.
Set up Sine-Gordon equation with initial and boundary conditions:
eq = {D[u[t, x], {t, 2}] == Sin[u[t, x]] + D[u[t, x], {x, 2}],
u[0, x] == E^(-x^2), Derivative[1, 0][u][0, x] == 0,
u[t, -10] == u[t, 10]};
Solve it:
sol = NDSolve[eq, u, {t, 0, 30}, {x, -10, 10}];
Plot it:
...
Only top voted, non community-wiki answers of a minimum length are eligible
