Hot answers tagged numerical-integration
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
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. ...
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
7
I think your initial condition is singular. In order to solve the ODEs, consider $(x(t), y(t))$ as a planar curve, we may try changing the parameter $t$ to the arc length parameter $s$:
$$\left\{\begin{split}
\frac{\mathrm{d}x}{\mathrm{d}s}=\frac{x'(t)}{\sqrt{x'(t)^2+y'(t)^2}}\\
\frac{\mathrm{d}y}{\mathrm{d}s}=\frac{y'(t)}{\sqrt{x'(t)^2+y'(t)^2}}\\
...
6
This will work:
\[Beta] = \[Omega]0/4; \[Omega] = 2 \[Pi]; \[Omega]0 = 3/2 \[Omega];
tend = 500;
s = ParametricNDSolveValue[{x''[t] +
2 \[Beta] x'[t] + \[Omega]0^2 Sin[
x[t]] == \[Gamma] \[Omega]0^2 Cos[\[Omega] t], x[0] == 0,
x'[0] == 0}, x, {t, 0, tend}, \[Gamma], MaxSteps -> Infinity];
Manipulate[
Plot[Evaluate[D[s[\[Gamma]][t], ...
6
This is a classical shock-tube problem in which a initially diaphragm separates a hi-pressure, high-density region from one of lower pressure and density. The classical exact solution has multiple discontinuities, a shock wave and a contact-surface (density discontinuity) that propagate to the right, and a continuous rarefaction wave traveling into the ...
6
As has been noted by ruebenko in the comments, there does seem to be a bug in the handling of infinite-range Bessel function integrals when MinRecursion and MaxRecursion are both set to non-default values. For instance, even the simple
NIntegrate[BesselJ[0, x], {x, 0, ∞}, MinRecursion -> 10, MaxRecursion -> 15]
chokes with a NIntegrate::minmax ...
5
The problem is that Mathematica prematurely threads r[t] - p not knowing r[t] is actually in $\mathbb{R}^2$
In[]:= r[t]-{0,0}
Out[]= {r[t],r[t]}
Which is not what you want. A quick fix for these types of issues is to create a function that only evaluates for numerical values (Changed to NDSolve since I only have v8):
dummy[r_?(VectorQ[#, NumericQ] ...
5
- update -
@whuber gives insightful comment and I agree I should mention here the following.
Below we consider a perturbed form (solvable exactly in Bessel functions) of original equation up to linear term. This allows us to understand behavior of the system around t ~ 0.
- original -
The following shows that under some general assumptions there is no ...
4
This integral equation is solvable using the LaplaceTransform technique:
Clear[s, t];
eqn = y'[t] == -Integrate[y[t1] Exp[t1 - t], {t1, 0, t}]
LaplaceTransform[eqn, t, s]
(*
==>
s LaplaceTransform[y[t], t, s] - y[0] == -(
LaplaceTransform[y[t], t, s]/(1 + s))
*)
Solve[%, LaplaceTransform[y[t], t, s]]
(*
==> {{LaplaceTransform[y[t], t, s] -> ...
4
The angle signs mean that you're looking at a shorthanded expression. Take a look at
FullForm[s]
to see the real thing.
The curly braces are there because NDSolve is going to give you consistently formatted results whether you look for only one function with only one solution or your solution involves more than one function, and more than one possible ...
3
You can define a solution function depending on your parameter and then use it for the animation :
sol[sigma_] := {x, y, z} /.
NDSolve[{x'[t] == sigma (y[t] - x[t]), y'[t] == 28 x[t] - y[t] - x[t] z[t],
z'[t] == x[t]*y[t] - 8/3 z[t], x[0] == z[0] == 0, y[0] == 2},
{x, y, z}, {t, 0, 25}][[1]]
Animate[
With[{f = sol[sigma]},
...
3
An alternative approach is to form an approximate interpolating function from your actual function, which is (often) cheaper to integrate. To wit,
nsin = FunctionInterpolation[Sin[x], {x, 0, 10}];
ni = Derivative[-1][nsin];
Plot[ni[t], {t, 0, 10}]
In this case, we know that the integral of sine can be expressed simply, so we can compare the exact ...
3
Using exact values (e.g. setting γ = 1/5 instead of γ = 0.2) and increasing the MaxSteps option should yield a reasonable result. Below, I add the filigree of also using a method I tend to prefer for problems with oscillatory solutions (although I've found that just setting Method -> "StiffnessSwitching" without increasing MaxSteps still works well):
xp ...
3
By default, NDSolve uses a "shooting" method to satisfy the boundary conditions: it picks an initial condition and then evolves the equation to see what boundary values are produced by the initial conditions. It then rejiggers these initial conditions to produce the boundary values specified by the boundary conditions. Here is the Mathematica documentation ...
3
"LevinRule" should work splendidly here, I think:
NIntegrate[-m Exp[-m] BesselJ[1, m]^2, {m, 0, Infinity},
Method -> "LevinRule", WorkingPrecision -> 20]
-0.18196415067209554877
ruebenko's answer has given a closed form for this particular definite integral. Personally, I prefer it when the parameters of the elliptic integrals are ...
3
@ruebenko confirmed that this is a bug.
Here's a workaround:
currentTimeStep = 0;
monitor[t_?NumericQ] := (currentTimeStep = t; 1)
pfun = ParametricNDSolveValue[{y'[t] == a monitor[t] y[t], y[0] == 1},
y, {t, 0, 10000}, {a}, MaxSteps -> 100000];
Dynamic[currentTimeStep]
pfun[1]
2
This is too long for a comment, but it isn't an answer.
Perhaps you would like to consider a more compact way to write down your equations:
m = -(6.672*10^-11) (7*10^17) ;
st = {{1000, 1000, 1000, 0, -100, 0},
{500, -1000, -1000, -110, 100, 0},
{0, 100, 500, 350, -100, 0}};
r = {x @ t, y @ t, z @ t};
o[n_] := NDSolve[Join[{Equal @@ ...
2
Your functions for y and x are single values and not lists, and thus doing a ParametricPlot for the sum of these is not going to work. However, to solve the equations and get a plot you can run the following:
alpha = 0.5
beta = 1
solf[b_, c_, tp_] := Module[{sol},
sol = NDSolve[{x'[t] == -x[t]*a[t] - 3*x[t] + alpha*y[t]^2 +
1/2*beta*(r1[t]^2 - ...
2
As an illustration of my initial comment, let's look at various starting points, following @Andrew's method:
sols = Table[
NDSolveValue[(2 + Sqrt[2] + s^2 - (-2 + Sqrt[2]) s^4 + s^6)/(1 + s^2)^2 yy[s] -
(1 - s^2)^2 yy''[s] == 0 && yy[start] == yy'[start] (start + 1) && yy'[start] == 1,
yy[s], {s, start, 1}], {start, -1 + ...
1
As a workaround it is possible to step up a bit from the left end:
leftend = -1 + 10^-8;
NDSolve[(2 + Sqrt[2] + s^2 - (-2 + Sqrt[2]) s^4 + s^6)/(1 + s^2)^2 yy[s]
- (1 - s^2)^2 yy''[s] == 0 && yy[leftend] == 0 && yy'[leftend] == 1, yy[s],
{s, leftend, 1}][[1]]
1
Here is a way:
sol = ParametricNDSolveValue[{-y*a'[x] == (a[x]^2 - 1)*x + a[x],
a[-10] == 0.6}, a, {x, -10, 10}, {y}]
sol[1]
NIntegrate[sol[y][t] /. y -> 1, {t, 0, 1}]
or alternatively via a Table:
ListLinePlot[
Table[NIntegrate[sol[y][t], {t, 0, 1}], {y, 0.1, 1, 0.1}]]
1
The issue is not with NDSolve but plotting:
\[Beta]=\[Omega]0/4;\[Gamma]=1/5;\[Omega]=2 \[Pi];\[Omega]0=3/2 \[Omega];
tend = 500;
s = NDSolveValue[{x''[t] +
2 \[Beta] x'[t] + \[Omega]0^2 Sin[
x[t]] == \[Gamma] \[Omega]0^2 Cos[\[Omega] t], x[0] == 0,
x'[0] == 0}, x, {t, 0, tend}, MaxSteps -> Infinity]
Plotting a smaller subrange helps:
...
1
If you just want to do a simple cartoon (as opposed to evaluating the solution components of the Lorenz equations at particular values), you can just directly extract the points generated by NDSolve[]. Here's one way to go about it:
sol[σ_] :=
Transpose[Through[{x, y, z}["ValuesOnGrid"] /.
First @ NDSolve[{x'[t] == σ (y[t] - x[t]),
...
Only top voted, non community-wiki answers of a minimum length are eligible
