Tag Info

Hot answers tagged

25

Generally always check Demonstrations site for good code. I cannot not mention an excellent "classic" of planar three body problem by Stephen Wolfram and Michael Trott. Code is short and I verified it runs in the latest M10.1. I slightly changed variable labels so code parses better here, removed MaxRecursion -> ControlActive[3, 9] from plot option and ...


8

eqns = { x1'[t] == -10*x1[t] + 10*x2[t], x2'[t] == a1*x1[t] - x2[t] - x1[t]*x3[t], x3'[t] == -8*x3[t]/3 + x1[t]*x2[t], x1[0] == 1, x2[0] == 1, x3[0] == 1}; Plot[ Evaluate[x1[t] /. Table[ NDSolve[eqns, {x1[t], x2[t], x3[t]}, {t, 0, 10}, MaxStepFraction -> 1/10000][[1]], {a1, {50, 55}}]], {t, 0, 10}, PlotStyle -> ...


8

Here's a slightly different approach, using ParametricNDSolve: sols = ParametricNDSolve[ {x1'[t] == -10 x1[t] + 10 x2[t], x2'[t] == a1 x1[t] - x2[t] - x1[t] x3[t], x3'[t] == -(8/3) x3[t] + x1[t] x2[t], x1[0] == 1, x2[0] == 1, x3[0] == 1}, {x1, x2, x3}, {t, 0, 10}, {a1} ]; This produces ParametricFunction objects, which take a value ...


7

The reason why you get a factor of 500 of I have explained in my comment. Let's replace your g with a better behaved function: mass = 100; width = 10^-2; g[x_] := mass/width HeavisideLambda[x/width] This is a triangular peak with area 100 and basewidth 0.01. Now let's impose the desired boundary conditions of $\partial u / \partial x =0$ at the ...


7

Plot[Evaluate[Through@pfun[50, 1, 1, 1]@t], {t, 0, 10}] (* or Plot[#[t] & /@ pfun[50, 1, 1, 1], {t, 0, 10}, Evaluated -> True] *) Note: The parametric function pfun[50, 1, 1, 1] is a list of 3 InterpolatingFunctions: pfun[50, 1, 1, 1] {InterpolatingFunction[{{0., 10.}}, <>], InterpolatingFunction[{{0., 10.}}, <>], ...


6

In a sense, the returned solutions are valid, but only over a restricted domain. We can determine the domain for each solution with something like Reduce[ode && -Infinity < t < Infinity /. sol, t] We can join the domain with the corresponding solution via ConditionExpression. ode = {x'[t] == Sqrt[x[t]], x[0] == 4}; dsol = DSolve[ode, x, ...


5

Correct sin/cos to Sin/Cos. Define your functions properly: w1[y_] :=.... Then put your conditions to Solve and done. w1[y_] := (E^(β*y))*(A1*Cos[β*y] + A2*Sin[β*y]) + (E^(-β*y))*(A3*Cos[β*y] + A4*Sin[β*y]); w2[y_] := (E^(-β*y))*(A5*Cos[β*y] + A6*Sin[β*y]); vars = Symbol["A" <> ToString[#]] & /@ Range[6] {A1, A2, A3, A4, A5, A6} ...


5

First, you're not using a fixed step method. (An Euler scheme may be applied to any step size and to one that varies.) To get a true fixed step method you have to turn off "DiscontinuityProcessing" when you have a discontinuous ODE; otherwise, NDSolve will try to adapt the steps to account for the discontinuity. The "DiscontinuityProcessing" stage resets ...


4

Update First I'd like to mention that after checking the definition of bilateral Laplace transform and Fourier transform carefully, I'm sure currently the formula for the relationship between them on the wikipedia page is wrong, the correct one should be: $$\mathcal{B}\{f(t)\}(s) = \sqrt{2\pi}\mathcal{F}\{f(t)\}(is)$$ Then we can use the internal function ...


4

The reason it's not finishing is that you set MaxSteps -> Infinity with a high AccuracyGoal. We can actually solve this system analytically, by replacing NDSolve with DSolve: qq = DSolve[{q'[t] == v[t]/r - q[t]/r*(1/c), q1'[t] == (((q[t] + q1[t])/(2*a*ϵ0*k2))*δ), q[0] == 10*10^-9, q1[0] == 10*10^-9}, {q, q1}, t]; We can inspect the solution: ...


4

By request of OP DE = p'[h] == -k p[h]; (* Assign an equation to the symbol DE *) (* semicolon omits the output *) dsol = DSolve[{DE, p[0] == 1013}, p[h], h] (* Assign the rule that is produced for a p[h], such that it would be the solution to ...


4

[Edit notice: It turns out the system can be integrated exactly, which leads to faster performance.] Using DSolve I replaced the approximate coefficients 0.2 and 0.8 by the exact numbers 2/10 and 8/10; otherwise, things that should cancel out do not and lead to false singularities. myode1 = Function[{kp1, kp2, θ, temp, HRT, p1in}, Evaluate@ ...


4

As noted in my comment, {x, y, z} must be defined in order for ItoProcess to work. Even then, the solution is unstable for the parameters chosen in the question. Stable parameters are, for instance, sigma = .01 and a RandomFunction step size of 0.1. In all, the modified code sigma = .01; trajectory = ItoProcess[{\[DifferentialD]v[t] == (a*δ*to[t] - ...


4

This can actually be solved analytically (DSolve), noting it is only algebraic in t: urt[r_, t_] = Simplify[(u /. First@DSolve[{kh t D[u[r], r] - 2 Pi a u[r] + 2 Pi a == 0, u[302] == 21/100}, u, r])[r]] 1 - 79/100 E^((2 a Pi (-302 + r))/(kh t)) Plot3D[urt[r, t] /. {kh -> 10^-6, a -> 300}, {t, 10, 200}, {r, 298, ...


3

This is a non linear equation which does not seem to have a closed form solution. Try instead something like this: eqn1 = D[y[t], {t, 2}] - (-M Cos[y[t]] Sin[y[t]] (y'[t])^2 - g/l Sin[y[t]])/(1 - M (Cos[y[t]])^2); eqn2 = D[x[t], {t, 2}] - (M g Sin[y[t]] Cos[y[t]] + M l Sin[y[t]] (y'[t])^2)/(1 - M (Cos[y[t]])^2); eqnSet = {eqn1 == 0, eqn2 ...


3

The differential equation, its initial condition, and its boundary conditions are translationally invariant in space. Consequently, the solution must be independent of x and y. Indeed, solving the equations as given in the Question does give a spatially constant solution that oscillates in time. For instance, DensityPlot[Evaluate[Re[A[x, y, 10000]] /. ...


3

I just wanted to point out that all the routes above will work if there is only one differential equation: data = NDSolveValue[{ x''[t] - k1*(1 - x[t]^2)*x'[t] + k2*x[t] == 0, x[0] == 2, x'[0] == 0} /. {k1 -> 1., k2 -> 1.}, Table[{t, x[t] + RandomReal[{-.3, .3}]}, {t, 0, 10, .2}], {t, 10}]; dataT = data\[Transpose]; ti = dataT[[1, ...


3

There is a problem with your initial velocity: Hmax = Vo^2*Sin[phi Degree]^2/(2*9.8) (*13.9274*) With that initial Velocity, you can only reach 13.9 meters (far from the 25+ from your data), and that's without drag. So you can't fit anything really. Maybe the velocity is not well measured. What you could do is fit Cdrag and Vo. Using part of Mariusz's ...


3

The FullForm of is Derivative[1][y][x] so when Mma goes to match the pattern y[x], there is no match. I'll use the code-simplifying assumption that you want y[x] = x^2. The issue is that y'[x] /. y[x]->x^2 evaluates to Derivative[1][y][x]. Instead, try y'[x] /. y->Function[x, x^2] which evaluates to 2 x as expected.


3

Why not {sol} = NDSolve[{f'[x] == g[x], g'[x] == f[x], f[0] == 1, g[0] == 0}, {f, g}, {x, 0, 1}] omitting the definition of f? Check: {Inactivate@NIntegrate[g[x0], {x0, 0, x}], f[x]} /. sol /. x -> 0.5 // Activate (* {0.521095, 0.521095} *)


2

You can use ParametricNDSolveValue for this task in two different ways. First, and more convenient way, is to use f''[val] as the second argument of ParametricNDSolveValue to get a function of the parameter a that gives the value of the function f'' at t==val {pf0, pf1, pf5} = ParametricNDSolveValue[{f'''[x] + f[x]*f''[x] == 0, f[0] == a, ...


2

I have changed so solution and vector plot match: Subscript[\[Tau], yz] = -\!\( \*SubscriptBox[\(\[PartialD]\), \(x\)]\(\[CapitalPhi][x, y]\)\); Subscript[\[Tau], xz] = \!\( \*SubscriptBox[\(\[PartialD]\), \(y\)]\(\[CapitalPhi][x, y]\)\); vp = VectorPlot[{Subscript[\[Tau], xz], Subscript[\[Tau], yz]}, {x, 0, 1}, {y, 0, 1}]; \[CapitalOmega] = ...


2

Starting from the ODE: DE = p'[h] == -k p[h]; Let us solve for the function h sol = First[DSolve[{DE, p[0] == 1013}, p, h]] Note that sol is a rule defining the function p; Then applying the corresponding rule, we can solve for k: rule2= First[Solve[ReplaceAll[p[20], sol] == 50 , k, Reals]]; ReplaceAll[ k, rule2] Now we can define a new ...


2

What follows is my latest edit. I have moved my previous attempts to solve this to the bottom of the answer. EDIT 04.05.15 I have given this problem some more thought. It's peculiar, that Mathematica cannot solve it straight away, but it should in fact be quite easy. Using your notation, we are dealing with the following differential equation: $$ ...


2

Define this function: f[k1_?NumericQ, k2_?NumericQ, k3_?NumericQ] := Sum[Total[(ci[[i, All]] - Map[pfun[k1, k2, k3][[i]], ti])^2], {i, 1, 3}] // Quiet Then, fit = NMinimize[f[k1, k2, k3], {k1, k2, k3}]; params = fit // Last (*{k1 -> 0.000194805, k2 -> 0.0291469, k3 -> 0.109229}*) Plot it, Table[Show[ ListPlot[Transpose[{ti, ci[[i]]}]], ...


2

ClearAll[Cto, kd1, ke2, kf3, FA, FB, FC, FD, FE, FF, Fto, V] Cto = 3.5; kd1 = 0.25; ke2 = 0.1; kf3 = 5.0; Note: = (not ==) in the line above. s = NDSolve[{FA'[V] == -kd1*Cto^3*(FA[V]/Fto[V])*(FB[V]/Fto[V])^2 - 3 ke2*Cto^2*(FA[V]/Fto[V])*(FD[V]/Fto[V]), FB'[V] == -2 kd1*Cto^3*(FA[V]/Fto[V])*(FB[V]/Fto[V])^2 - ...


2

Here is one variant that seems to evaluate to a number. myode1[p1in_?NumericQ, HRT_?NumericQ, temp_?NumericQ] := ParametricNDSolveValue[{p1'[t] == 1/HRT (0.2 p1in - p1[t]) - kp1*Exp[θ (temp - 20)] p1[t], p2'[t] == 1/HRT (0.8 p1in - p2[t]) - kp2*Exp[θ (temp - 20)] p2[t] + kp1*Exp[θ (temp - 20)] p1[t], p1[0] == 0.2 p1in, p2[0] ...


2

First, with expressions like these E^(-(0.00121935/(q[t] + q1[t]))) ... E^(-((0.00121935 Sqrt[1 + 24213.1 (q[t] + q1[t])]) it's unlikely that the ODE can be integrated exactly, as has been pointed out in comments. Second, using approximate coefficients (instead of exact rational, algebraic, or transcendental numbers) is going to make rather difficult ...


2

If I change one of your coefficients, I get the example plot: equation1 = {x'[ t] == (1/(1 + (5*y[t] - 5*(2/(1 + Sqrt[1 + 1600*y[t]]))*y[t])^2)) - 0.1*x[t], y'[t] ==(*0.05*)0.5*x[t] - (* <-- N.B. *) ((10*y[t]*(2/(1 + Sqrt[1 + 1600*y[t]])) + 0.03*y[t])/(0.05 + y[t])) - 0.1*y[t], ...


2

There is a DSolve and DSolveValue function for this. In case of $m$, $b$ and $k$ constant it can be easily solved: DSolveValue[m x''[t] + b x'[t] + k x[t] == DiracDelta[t], x[t], t] (* E^(1/2 (-(b/m) - Sqrt[b^2 - 4 k m]/m) t) C[1] + E^(1/2 (-(b/m) + Sqrt[b^2 - 4 k m]/m) t) C[2] - ((E^(1/2 (-(b/m) - Sqrt[b^2 - 4 k m]/m) t) - E^( 1/2 (-(b/m) ...



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