New answers tagged

1

Oh! I'm coming to late! r1 = 1.8; r2 = 2.2; ode1 = h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r]); ode2 = p'[ r] == -(2 p[r]/r) - (h'[r]/h[r])* p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r])); sol = ParametricNDSolve[{ode1, ode2, m[r] == (r/2)*(1 - h[r]), m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}] One way to find ...


1

The second comment on this post is incorrect. You can set a condition for any time t0, whether initial or not. Independently of this, you can solve the equation for any interval of the variable. Example: NDSolve[{f'[t] == -f[t], f[2] == 1}, f, {t, 1, 3}] Here the equation is solved for the interval $t \in [1,3]$, but the condition is set for $t == 2$ as ...


2

r1 = 1.8; r2 = 2.2; ode1 = {h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r])}; ode2 = {p'[r] == -(2 p[r]/r) - (h'[r]/h[r])* p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]))}; hpm = ParametricNDSolveValue[{ode1, ode2, m[r] == (r/2)*(1 - h[r]), m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}] Plot[Evaluate[Through[hpm[0.8, 1.][t]]], ...


5

It will work perfectly once you write it correctly x /. NSolve[{3 BesselJ[1, x] + x D[BesselJ[1, x], x] == 0, 16 > x >= 0}, x] {2.9496, 5.84113, 8.87273, 11.9561, 15.0624} You made a mistake when you write D[BesselJ[1, x]] and also you are missing 3 in first term.


0

You have your NDSolve xsoln = NDSolve[{x'[ t] == -Sqrt[(1/(4 - 2 Sqrt[2]) + (4 - 2 Sqrt[2])/(2 x[t]))^2 - 1], x[0] == 4}, x, {t, 0, 10}] which returns a replacement rule. This replacement rule will work with derivatives, just try Plot[{x[t], x'[t]} /. xsoln, {t, 0, 10}, Evaluated -> True] to see this. Then you solve the algebraic ...


4

Here is a simple example for the default method (LSODA) that shows a couple of the issues related to the situation in the question: $$y' = \exp(-10 \, x^2), \quad y(-20) = 1, \quad -20 \le x \le 20$$ The exact solution is a scaled error function (Erf), which has a sigmoid graph going from $y(-20)=1$ to $y(20) \approx 1.5605$. For $|x|$ large, $y'$ is ...


8

How can I constrain the locator to stay within the region defined by RegionPlot? You can check if locator's coordinates fulfill the condition defining your region. It can be done with the second argument of Dynamic if you introduce Locator explicitly. Take a look at line with Locator[Dynamic[p, With[{... A strange behavior occurs whenever the left ...


2

Here is what I get when I try your equation: a = 1; b = 0.2; c = 4; y0 = 8.5; v0 = 0; DSolve[{a y''[t] + b y'[t] + c y[t] == 0, y[0] == y0, y'[0] == v0}, y[t], t] (* {{y[t] -> 8.5 E^(-0.1 t) (1. Cos[1.9975 t] + 0.0500626 Sin[1.9975 t])}} *) You may need to restart your kernel and try again.


5

x = FindInstance[{Cot[Sqrt[λ]] == Sqrt[λ], λ >= 0, λ < 100}, λ, 100] You can take the numerical value of the solution using N and use ReplaceAll (/.) to get the values and finally take the square root: λ /. N[x] // Sqrt {0.860334, 3.42562, 6.4373, 9.52933} This should work regardless of the number of solutions.


2

I'm on MMA 10.3 on OSX 10.10.5 Running your code indeed gives me the same result, a straight line. However this is not surprising to me as Evaluate[y /. nds] returns, unexpectedly, y Therefore in the plot we are plotting, y from 0 to 30, hence a straight line, gradient 1, from 0 to 30... Furthermore the syntax in your NDSolve is a little off to ...


12

Some explanations first The substitution in the question introduces the reduced wave function $u(r)$ by solving the original radial equation in polar coordinates, $$-\frac{1}{2}\left(R''(r)+\frac {1}{r}R'(r)\right) - \frac{1}{r}R(r) + \frac {m^2}{2r^2}R(r) = E R(r)$$ using the ansatz $$R(r)\equiv \frac{1}{\sqrt{r}}u(r)$$ The apparently divergent ...


3

Below is a workaround for the simple case. The OP can say whether it works more general. I haven't quite tracked down yet why the system is set up incorrectly with the default Method -> {"EquationSimplification" -> "Solve"} and with Method -> {"EquationSimplification" -> "Residual"}. But it works in this case with Method -> ...


2

happy fish ,he said "nth derivative is not natively supported" Yes it's true,but from here. Method1: For simple functions you can use InverseFourierTransform. f[s_] := Sin[s]; nthDeriv1[f_, s_, n_] := FullSimplify[InverseFourierTransform[(-I k)^n FourierTransform[f, s, k], k, s], {n \[Element] Integers, n > 0}] nthDeriv1[f[s], s, n] $$\sin ...


0

A method to solve the original decoupled system : s = NDSolve[{ D[u[x, t], t] == -D[u[x, t], x] - u[x, t], u[x, 0] == 1, u[0, t] == 1}, {u[x, t]}, {x, 0, 4}, {t, 0, 1}] v[x_, t_] = (u[x, t] /. s[[1, 1]]) /. x -> x + 1; w[x_, t_] = (u[x, t] /. s[[1, 1]]) /. x -> x + 2; z[x_, t_] = (u[x, t] /. s[[1, 1]]) /. x -> x + 3; Row[{ ...


1

As suggested by the Documentation Center, you can remove excess (last 3) conditions. s = NDSolve[{D[u[x, t], t] == -D[u[x, t], x] - u[x, t], D[v[x, t], t] == -D[v[x, t], x] - v[x, t], D[w[x, t], t] == -D[w[x, t], x] - w[x, t], D[z[x, t], t] == -D[z[x, t], x] - z[x, t], u[x, 0] == 1, v[x, 0] == 1, w[x, 0] == 1, z[x, 0] == 1, u[0, t] ...


4

Changing parameter values during integration works better with DiscreteVariables. But I think the problem with OP's code, in the question and the OP's answer, has more to do with Mathematica numerics. My solution Clear[bind]; zdot = 1/2 (1 - z[t]); ydot = 1/20*y[t] + z[t] - x[t]; xdotbind = D[Solve[-ydot - zdot == 0, x[t]][[1, 1, 2]], t] /. {y'[t] -> ...


4

If you look at the FullForm of your variable result you will see that it contains some small complex numbers. You can remove these by using Plot[Abs[result[[i]]] instead.


2

So, after playing around with this problem all day, I finally stumbled upon a version of the code that gives me my desired result, although I'm not at all clear on why it works where other versions fail. First, here is the working version: zdot=.5*(1-z[t]); ydot=.05*y[t]+z[t]-x[t]; ...


7

This might be useful if your actual event is more complicated, use a discrete variable as a flag: sols = NDSolve[{(1 + 25 Exp[-u[s]^2] u[s]^2) u'[s]^2 == u[s]^2, u[0] == 0.5, a[0] == 1, WhenEvent[u[s] == 4, a[s] -> 0], WhenEvent[a[s] == 0, "StopIntegration"] }, {u, a}, {s, 0, 20}, DiscreteVariables -> {a[s]}] u /. ...


4

There is no solution to your equation. If you specify Reals for k Mathematica outputs an empty list as solution: NSolve[113.68/k + 1345.6/k^2*(1 - Exp[-k/11.6]) == -4.9, k, Reals] (*{}*) If you plot the righ and left hand side of the equation they don't cut each other: Plot[{113.68/k + 1345.6/k^2*(1 - Exp[-k/11.6]), -4.9}, {k, -60, 60}]


4

The corrected equation can be solved using FindRoot: FindRoot[(113.68 t/k - (1345.6/(k^2))*(1 - Exp[-k*t/11.6]) == -4.9 t^2) /. t -> 1, {k, 0.5}] (* {k -> 0.235142} *) This solves the equation for the value of $k$ corresponding to $t$ = 1 second. FindRoot requires an initial "guess" for the value of $k$; in this case, I've used $k = 0.5$. Other ...


6

The problem seems to be that you give the initial condition at a different point that your lower limit of the integration. If you change x1 by 0 things seem to work: (*INPUT*)MP = 1; m = 0.2 MP; c = Sqrt[3/2]; x1 = 0; x2 = 4000; (*ODE*) Clear[s2] Table[s2[i] = NDSolve[{y''[x] + c (Sqrt[m^2 y[x]^2 + y'[x]^2]) y'[x] + m^2 y[x] == 0, y[0] == ...


2

Thanks to those that responded. Applying the boundary in the pde as a formula is okay for simple geometries but for more complicated ones it is cumbersome. I have tried user21's first suggestion: bmesh = ToBoundaryMesh[DiscretizeGraphics[heatsink], "BoundaryMarkerFunction" -> (2 & /@ # &)]; emesh = ToElementMesh[bmesh]; ...


4

I am not sure I understand the question 100% but here is what I think you are looking for: {vals3DL0, funs3DL0} = NDEigensystem[-(R''[r]/2) - (R'[r]/r) - (R[r]/r), R[r], {r, 0, 200}, 3, Method -> {"Eigensystem" -> {"Arnoldi", "Criteria" -> "RealPart"}, "SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" -> \ ...


2

A bit long for a comment: What I'd do is split all polygons into quad elements and use those directly. Roughly like this: sidesback1 = Partition[Table[i, {i, 1, nodes}][[3 ;; -1]], 4]; sidesfront1 = Partition[Table[i, {i, nodes + 1, 2*nodes}][[3 ;; -1]], 4]; bmesh = ToBoundaryMesh["Coordinates" -> Join[face1, face2], "BoundaryElements" -> { ...


3

soln = ParametricNDSolve[{x'[t] == -I a x[t] + I Sqrt[2] b y[t], y'[t] == I Sqrt[2] b x[t] + I Sqrt[2] b z[t], z'[t] == I Sqrt[2] b y[t] - I a z[t], x[0] == 1, y[0] == 0, z[0] == 0}, {x, y, z}, {t, 0, 20}, {a, b}] Plot[Through[{Re, Im}@#] & /@ {x[.1, .5][t], y[.1, .5][t], z[.1, .5][t]} /. soln, {t, 0, 10}, Evaluated -> True, ...


3

Equations use Equal (==) vice Set (=) eqns = {x'[t] == -I a x[t] + I Sqrt[2] b y[t], y'[t] == I Sqrt[2] b x[t] + I Sqrt[2] b z[t], z'[t] == I Sqrt[2] b y[t] - I a z[t], x[0] == 1, y[0] == z[0] == 0}; soln = DSolve[eqns, {x, y, z}, t][[1]] // FullSimplify {x -> Function[{t}, (4*b^2*E^((-I)*a*t - (1/2)*I* (a + Sqrt[a^2 ...


2

There are two different methods for solving the OP"s problem : The Method of Lines with the option "SpatialDiscretization" -> {"TensorProductGrid"... The Method of Lines with the option "SpatialDiscretization" -> {"FiniteElement"}. This solution is the Mathematica 10 implementation of the Finite Element Method for transcient PDEs. In both cases ...


2

A point has measure 0 in any numeric approximation possible, so it's not only mathematica who doesn't like point-like conditions. soln = NDSolveValue[{Laplacian[u[x, y], {x, y}] == NeumannValue[-1, y == -1. && Abs[x] < 0.1] + NeumannValue[1, y == 1 && Abs[x] < 0.1], u[0, y] == 0}, u, {x, y} \[Element] Rectangle[{-1., ...


2

The built in functions Solve, NSolve, DSolve, etc. return a list of solutions, even if there is only one. If you are sure that you want only a single solution you could use y1[x_] = First[y[x] /. DSolve[{y''[x] + 0.25 y'[x] + y[x] == 0, y[0] == 1/2, y'[0] == 5/4}, y[x], x]] Then evaluating y1[x] will yield a single value The plot below was produced by ...


2

Grad (g) returns a vector - but your PDE is scalar. Here is an example where I sum over the squares g and take the sqrt (which is a scalar): Needs["NDSolve`FEM`"]; reg = ImplicitRegion[x^2 + y^2 <= 1, {x, y}]; mesh = ToElementMesh[reg]; (*RegionPlot[reg]*) s = NDSolveValue[{Derivative[0, 2][f][x, y] + Derivative[2, 0][f][x, y] == 0, ...


3

I would memoize sol[a]. The Evaluate in position does nothing if it does not wrap the entire expression after the :=. It's not that important, so I would just drop it. The issue with [[1]] (or First) can be handle in sol. Here are the changes I've described: sol[a_] := sol[a] = First@NDSolve[{x''[t] == -2 x[t], x[0] == a, x'[0] == 1}, {x}, {t, 0, 10}] ...


0

This helps to get rid of that [[1]] as per your request, however without having the code that is really slowing this down, I cannot be sure it would help improve your performance. If it does, great! If not, you have eliminated one possible cause. You state the performance glitch is with position[t,a] function. However, can you remove the 0.01 timesteps in ...


3

I'm answering my own question since it may eventually help somebody else who is facing a similar problem. In this case since the equations are linear on the particles' velocities one can find "easily" the one that makes the force zero. This is how I'm doing it. First, let's put the problem in abstract terms and for only one dimension for the sake of ...


2

One could also feed the ODE into DifferentialRoot[] (along with arbitrary initial conditions) and then apply FunctionExpand[]. sol = FunctionExpand[DifferentialRoot[Function[{y, x}, {y''[x] + 2 y'[x] + 5 y[x] == x Cos[x], y[0] == C[1], y'[0] == C[2]}]][x]]; ...


2

On my mac doesn't work either. But you can do << NDSolve`FEM` l = 10.; w = 10.; h = 5.; reg = Cuboid[{-l/2, -w/2, -h/2}, -{-l/2, -w/2, -h/2}]; mesh = ToElementMesh[reg, MaxCellMeasure -> {"Area" -> .1}] (*ElementMesh[{{-5., 5.}, {-5., 5.}, {-2.5, 2.5}}, {HexahedronElement[ "<" 16384 ">"]}]*) Somehow ImplicitRegion works like crap ...


2

With version 10.4.1 on Linux I get << NDSolve`FEM` l = 10.; w = 10.; h = 5.; reg = ImplicitRegion[-l/2. <= x <= l/2. && -w/2. <= y <= w/2. && -h/2. <= z <= h/2., {x, y, z}]; mesh = ToElementMesh[reg, MaxCellMeasure -> {"Area" -> .1}]; mesh["Bounds"] {{-5.`, 5.`}, {-5.`, 5.`}, {-2.5`, 2.5`}} Other ...


10

Here is how you would do it using the standard add-on package VariationalMethods, which is meant for calculations like this: Clear[m, k, c, x, t]; T = 1/2 m x'[t]^2; V = 1/2 k x[t]^2; L = T - V; Needs["VariationalMethods`"] hamiltonianEq = h == FirstIntegral[t] /. Last[FirstIntegrals[L, {x[t]}, t]] (* ==> h == 1/2 (k x[t]^2 + m ...


5

genCoords = {x[t]}; ke = 1/2 m x'[t]^2; v = 1/2 k x[t]^2; q = -c x'[t]; l = ke - v; Solve for x'[t] in terms of p[t]: rule = First@Solve[p[t] == D[l, x'[t]], x'[t]] (* {x'[t] -> p[t]/m} *) and then replace this expression into the Hamiltonian: x'[t] D[l, x'[t]] - l /. rule (* p[t]^2/(2 m) + 1/2 k x[t]^2 *) Note that I have used lower-case symbols ...


4

Without diving into your code too much, everything will run a LOT more smoothly if you use ParametricNDSolve to solve your differential equations with the parameter a: pfun = ParametricNDSolveValue[{x''[t] == -2 x[t], x[0] == a, x'[0] == 1}, x, {t, 0, 10}, {a}] position[t_, a_] := {Sin[#], Cos[#]} &@pfun[a][t] You can keep everything else the same. ...


1

NDEigensystem was added in version 10.2 (or was it 10.3?) but version 10.0 is not going to work unless you use this.


2

As I noted in a comment above, it is quite possible for the solution of a nonlinear PDE to become singular at finite t, and that appears to be occurring here. And, as noted by Dr. Wolfgang Hintze, the right side of the PDE, when integrated over {x, -6, 6} is zero. So, the integral of f must be a constant, and indeed it is. For the parameters given in the ...


10

EDIT #2 My error was useful. It brought me to the conclusion that the difficulties in solving the PDE of the OP are due to the drift term $$\frac{\partial (x u(x,t))}{\partial x}$$ If the drift term is included, many boundary problems are ill defined. It turns out that there are cases where mathematically there is only a trivial solution u = 0 but ...


1

After some searching, it turns out that you can obtain those definitions using SubValues: Clear[f, x] f'[x] = 3 x SubValues[Derivative] (* Out: {HoldPattern[Derivative[1][f][x]] :> 3 x} *) As mentioned in this answer, SubValues is used for definitions of the following type: d[e][f] = something; As you know, f'[x] is really interpreted as ...


1

I don't know if it's a bug, but Mathematica makes a move that, when it's your student you say, why did you do that! It was simpler before. I ran this code without the replacement with Exp[stuff__ * C[1]] in order to figure out what to do to fix it. (The colored output below are from the Print statements showing the various calls to Solve that happen inside ...


0

As explained in the comment. The sol yields {y[t]->...}. This argument will be evaluated in Plot but not with Point. If the aim is just to see bounce on vertical line you just have to take first part of sol (e.g First[sol]. If you want to plot y v time he is one of many ways: sol = y[t] /. First@NDSolve[{y''[t] == -9.81, y[0] == 5, y'[0] == 0, ...


1

I corrected errors to get you started, but the physical equations themselves may not be accurate. springalt[c_, a_, x_, rs_, vs_, tmax_] := NDSolve[{ϕ'[t] == ω[t], ω'[t] == -c (((Sqrt[(Cos[ϕ[t]] + a)^2 + (Sin[ϕ[t]])^2] - 3)/ Sqrt[(Cos[ϕ[t]] + a)^2 + (Sin[ϕ[t]])^2])*a* Sin[ϕ[t]] + (((\[Sqrt]((x - Cos[ϕ[t]])^2 + ...


3

Oh boy, what a question! This is very similar to some stuff I played a few weeks ago (Kerbal, what a game!). What follows solves (I think) the question you are asking. An approach that seemed to to help was to split the problem into to: before, and after the engine burn. I do this with the knowledge that the most efficient landing will comprise just a ...


1

Maybe like this: sol = DSolve[{x''[t] + x'[t]/100 + x[t] == Exp[I*t], x[0] == 1, x'[0] == 1}, x[t], t] (* {{x[t] -> (1/39999) E^(-t/200) ((39999 + 3999900 I) Cos[(Sqrt[39999] t)/200] - 3999900 I E^((1/200 + I) t) Cos[(Sqrt[39999] t)/200]^2 - (19799 - 100 I) Sqrt[39999] Sin[(Sqrt[39999] t)/200] - 3999900 I E^((1/200 ...


3

Your system is singular and cannot be fixed with Adams's method (see also J.M.'s comment). You can solve it with DSolve. nsol = NDSolveValue[{y'[x] == y[x]*y[x]*Exp[x] - 2*y[x], y[0] == 2}, y, {x, 0, 10}] At x == 0.6931469688260267`, step size is effectively zero; \ singularity or stiff system suspected. {x1, x2} = nsol["Domain"][[1]] (* {0., ...



Top 50 recent answers are included