# Tag Info

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 ...

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], ...

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] -> ...

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] = ...

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

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], ...

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

- 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 ...

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] ...

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 ...

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, ...

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]

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

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 ...

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 -> ...

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

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, ...

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 ...

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}]

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

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 ...

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 ...

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

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

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] + ...

Only top voted, non community-wiki answers of a minimum length are eligible