New answers tagged differential-equations
0
As Andrew stated above, changing the lower limit to something non-zero takes care of the problem.
The following code gives a speedy answer:
With[{Dif = 2300}, pde = D[u[t, r], t] == Dif/r*D[u[t, r]*r, r, r] -
1/(0.0945 E^(0.000212 + 0.4077 r))*(u[t, r] + 1.2004);
soln = NDSolve[{pde, u[0, r] == 0, Derivative[0, 1][u][t, 10^-8] == 0,
u[t, 500] == 0}, ...
3
Ok, your case is the one which can be easily helped by memoization. Basically, it will serve to compute the set of values in sol in reverse order, de facto making it just a convenient device to turn your algorithm into an iterative one without changing it much.
All you have to do is to change your definition of sol as:
sol[n_] := sol[n] = ...
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]]
6
This method of developing a truncated solution can be done as below. I illustrate with an example that DSolve does not seem much to like.
ode = x''[t] - t*x'[t] + Sin[t] == 0;
initconds = {x'[0] == 1, x[0] == 0};
We create a differentail operator to create this ode.
odeOperator = D[#, {t, 2}] - t*D[#, t] + Sin[t] &;
Now set up our Taylor series as ...
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] ...
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}}\\
...
8
I think what he means is the following (I could be wrong though). In Mathematica white noise is represented by WienerProcess "also known as Brownian motion, a continuous-time random walk, or integrated white Gaussian noise." ~ Documentation. So adopting (actually simplifying) this example we get the following:
\[ScriptCapitalP] =
...
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 ...
10
A little bit more. Still not fully diagnosed, but the problem isn't due to DSolve
... :
s1 = DSolve[{x'[t] == f*x[t] (1 - (x[t]/b)) - l x[t]}, x[t], t];
s2 = DSolve[{x'[t] == e*x[t] (1 - (x[t]/b)) - l x[t]}, x[t], t];
And the problem shows up when matching the initial condition:
Solve[(x[t] /. s2[[1]] /. t -> 0) == 4/10, C[1]]
(*
{{C[1] -> ...
6
The problem can be reduced to the DSolve expressions:
DSolve[{x'[t] == a*x[t]*(1 - (x[t]/b)) - l*x[t], x[0] == 0.4}, x[t], t]
DSolve[{x'[t] == h*x[t]*(1 - (x[t]/b)) - l*x[t], x[0] == 0.4}, x[t], t]
One can see that alphabetical order appears important:
With[{a = Symbol@#},
Shallow @ DSolve[{x'[t] == a*x[t]*(1 - (x[t]/b)) - l*x[t], x[0] == 0.4}, x[t], ...
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 - ...
0
Could be interpreted as a nonlinear least squares problem.
pts = {{70.26, 45.78}, {71.04, 46.32}, {37.23, 24.67}, {37.91,
28.78}};
ndsoln[{k1_, k2_, k3_}, x0_] := x[720] /.
NDSolve[{x'[t] == -k1 x[t]^2 + k2 x[t] y[t],
y'[t] == k1 x[t]^2 - k2 x[t] y[t] - k3 y[t], x[0] == x0,
y[0] == 0}, {x, y}, {t, 0, 720}][[1]]
ssfun[{k1_?NumericQ, ...
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],
...
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
If this really can be expressed linearly, then expressing in Matrix form is going to be the easiest thing and the best from a computational perspective. For example, with your definitions:
R6 = kf*(z[1] + mu*z[5]) - 2*z[6] + mu*kd^2*z[7] - kd^2*z[6];
R7 = kf*(z[1] + mu*z[6] - 2*z[7]) + kd^2*z[8] + kd^2*z[7];
R8 = kf*z[1] + (mu*z[7] - 2*z[8]) + kd^2*z[9] + ...
1
Does this do the job?
poly=(kz py + py + kr kx)
variables = Variables[poly];
reprule =DeleteCases[If[StringTake[ToString[#], 1] == "k", # ->
Subscript[StringTake[ToString[#], 1],
StringTake[ToString[#], {2, -1}]]] & /@ variables, Null]
poly/.reprule
This will take all the variables in an expression, work out if they are a combination of k ...
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]
0
One option is to reduce the PrecisionGoal and AccuraryGoal. Those are usually set quite high and perhaps this is a viable option here. Setting
, PrecisionGoal -> 4, AccuracyGoal -> 4
for the Ns=10 case runs in 12 seconds, where the original code runs in 20 seconds.
It were good if you removed the calls to N and replaced the real numbers with exact ...
0
Now I tried to change the interval as follows:
A0 = 2.0377638272727268`;
A1 = -7.105521894545453`;
A2 = 9.234000147272726`;
A3 = -5.302489919999999`;
A4 = 1.1362478399999998`;
h0 = 45.5;
\[Sigma]M = 0.00592251;
\[Lambda]1 = 1.0253896074561006`;
\[Lambda]2 = 1.3079437258774012`;
f1 = A1 + 2 A2 y[x] + 3 A3 y[x]^2 + 4 A4 y[x]^3;
b = h0^2/12 (5 A1 + 8 A2 y[x] ...
2
There's a g in your differential equation, I think you actually mean g[y, x], right? This is the root of your problem. By the way, why do you make y an argument for functions f1, b, g? (It doesn't hurt though…) In fact, I think it will be conciser not to use function definitions in this case.
Then, your "StartingInitialConditions" for Shooting isn't so ...
2
Well, this answer is quite incomplete because I can't fix the problem. I post this answer just to point out that the true reason for the changeless sol is the BCs (D[u[x, y, t], y] /. y -> 0) == mu1[x] and (D[u[x, y, t], y] /. y -> 1) == mu2[x] are largely ignored by NDSolve. If you try:
Clear["`*"]
T0 = 500; a = 0.002; phi[x_, y_] = 300;
(* I ...
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] -> ...
2
When you try to find the general solution with initial condition y[x0]==y0 you see what happens
sol = y[x] /.
First@DSolve[{y'[x] - x^2 (y[x] + 1) (y[x] - 2)^2 == 0, y[x0] == y0}, y[x], x]
As you see when y0==2 is directly inserted into the solution the denominator becomes zero. At this point you could use Limit to investigate further.
2
There are two things about your Manipulate that stand out as possible sources of trouble.
The use of Module to localize variables within a Manipulate isn't a good idea. It is better to introduce an invisible control to create a local variable. Also note that the variable plt isn't needed at all.
For complex computations such as you are carrying out, it is ...
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 ...
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 ...
1
It seems to work if you replace Infinity with a smaller number :
s = NDSolve[{Derivative[3][f][x] + 1/2 f[x] Derivative[2][f][x] == 0,
f[0] == 0, f'[0] == 0, f'[#] == 1}, f, {x, 0, 1}] & /@ Range[1, 50, 5];
Plot[Evaluate[f[eta] /. s], {eta, 0, 1}, PlotRange -> All,
PlotLegends -> (ToString[#] & /@ Range[1, 50, 5])]
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 ...
2
We start with the given differential polynomials.
a1 = A1*D[p1[x, y, z], {z, 2}] -
A2*(D[p1[x, y, z], {x, 2}] + D[p1[x, y, z], {y, 2}] +
D[q1[x, y, z], {x, 1}] + D[q2[x, y, z], {y, 1}]);
a2 = A9*(A2*
D[(D[q1[x, y, z], {x, 1}] + D[q2[x, y, z], {y, 1}]), {x, 1}] +
A2*(D[q1[x, y, z], {x, 2}] + D[q1[x, y, z], {y, 2}])) -
A4*(D[p1[x, ...
3
As the other answers have mentioned, you can confirm the solution without specifying values for C[1] and C[2]
lhs = q''[x] + 2 x/(x^2 - 1) q'[x] - 4*q[x]/(x^2 - 1);
sol = DSolve[lhs == 0, q, x][[1, 1]];
FullSimplify[lhs /. sol]
0
When C[2] is set to zero, FullSimplify is unable to find the right transformations to reduce the expression to zero, and ...
1
You can do :
sol[x_] = q[x] /.
First@DSolve[q''[x] + 2 x/(x^2 - 1) q'[x] - 4*q[x]/(x^2 - 1) == 0, q[x], x]
FullSimplify[Derivative[2][sol][x] + 2 x/(x^2 - 1) Derivative[1][sol][x] -
4*sol[x]/(x^2 - 1)]
(* 0 *)
0
You do not need WhenEvent here as stated. You can use it though but you have to treat k as a function of t, that is a discrete state variable k[t].
Manipulate[
DynamicModule[{
sol = NDSolveValue[{
y'[t] == k[t] y[t],
y[0] == 10, k[0] == k0,
WhenEvent[t > 10, k[t] -> k[t] + s],
WhenEvent[t < 20, k[t] -> k[t] - s]
...
0
I think you may not chose one parameter to 0. E.g. try c2=2 and c1=1:
FullSimplify[
D[LegendreP[1/2 (-1 + Sqrt[17]), x] +
2 LegendreQ[1/2 (-1 + Sqrt[17]), x], {x, 2}] +
2*x/(x^2 - 1)*
D[LegendreP[1/2 (-1 + Sqrt[17]), x] +
2 LegendreQ[1/2 (-1 + Sqrt[17]), x], x] -
4*(LegendreP[1/2 (-1 + Sqrt[17]), x] +
2 LegendreQ[1/2 (-1 + Sqrt[17]), x])/(x^2 - 1)]
...
2
I'd like to extend the solution offered by Michael E2:
psol = ParametricNDSolve[{A'[t] == k1*A[t], A[0] == 10,
WhenEvent[t < 10, k1 -> (k1 + s)],
WhenEvent[t < 20, k1 -> (k1 - s)]},
A, {t, 0, 100}, {k1 \[Element] Reals, s \[Element] Reals}];
(*Plot[Evaluate[A[0.1, 0.2][t] /. psol], {t, 0, 30}]*)
Note that even in the case the DE ...
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}]]
2
You don't need WhenEvent[] for this:
Manipulate[Plot[(A/. NDSolve[{A'[t] == (k1 + s HeavisidePi[1/10 (t - 15)]) A[t], A[0] == 10},
A, {t, 0, 100}] [[1]] )[x], {x, 0, 30}],
{k1, 0.0, 1.0}, {s, 0.0, 1.0}]
1
I don't think you can use WhenEvent to do what you want. The value of k1 is passed in the DE in NDSolve, not the symbol. WhenEvent has the attribute HoldAll, so that it deals with k1 and s as Symbols. Perhaps you could use ParametricNDSolve (see below).
Perhaps you want something like this?
kparam[t_?NumericQ, k1_, s_] := If[10 < t < 20, k1 + s, ...
4
No "brute-force" playing with NDSolve, we can get an idea of attraction basins with the StreamDensityPlot and StreamPoints option in it.
Let's find e few points of interest where the vector flow becomes zero.
E.g.
{x, y} /. {ToRules @ LogicalExpand @
Reduce[y == 0 && -9 Sin[x] - 1/5 y == 0 && -5 < x < 10, {x, y}]}
{{0, 0}, {-Pi, ...
15
Here's a bruteforce way to do it for the simple case when the attractor is a fixed point.
Find a fixed point
Pick initial values for ODE
Solve ODE, see if it gets close to fixed point
Go back to 2 until satisfied
By looking at Reduce[y == 0 && -9 Sin[x] - 2/10 y == 0, {x, y}, Reals] we see that for instance {x=0,y=0} is a fixed point, let's use ...
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], ...
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:
...
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
For an n that "ensures x1 > x2 at all times", I think it may not be possible. Look at the plot below, it seems for t negative enough, x1 will eventually below x2 for all $1 \leq n \leq 8$:
Plot3D[
Evaluate[x1 - x2], {t, -4, 1}, {n, 1, 8},
PlotPoints -> 100,
MeshFunctions -> {#1 &, #2 &}, Mesh -> {{0}, 20},
MeshStyle -> ...
1
The problem is, that you call NDSolve while you haven't specified numerical values. Later, in your Table they are put in. In your second example you never say which value $\mu$ should get.
An easy way to fix this, is to make a function call out of your NDSolve. In this way it gets only evaluated when you put values in.
sol = Function[{x0, y0, mu},
...
Top 50 recent answers are included
