# Tag Info

14

You need a DirichletCondition (new in V10) here. Using regions (also from V10): Ω = RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]]; sol = NDSolveValue[{ D[u[x, y], x, x] + D[u[x, y], y, y] == 0, DirichletCondition[u[x, y] == 100., x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= ...

10

Here are a couple of suggestions. First let's look at the mesh generation with is documented with ToElementMesh and in the mesh generation tutorial. Needs["NDSolveFEM"]; x2 = 2; y2 = 1; r = 1/8; reg = ImplicitRegion[ 0 <= x <= x2 && 0 <= y <= y2 && (x - x2/2)^2 + (y - y2/2)^2 >= r^2, {x, y}]; For the mesh ...

9

The problem - Genericity of soutions The problem we encounter here is that DSolve can return only a general solution however that general solution cannot satisfy such an initial condition as y'[0] == 1. The issue is related to an arbitrary choice of constants of integration i.e. such constants that are specific to certain types of a differential equations ...

9

Here is what I think the issue is: Let's look at what NDSolve parses. Needs["NDSolveFEM"] {state} = NDSolve`ProcessEquations[{op == 0, bc}, u, {x, y} \[Element] reg, Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 0.005}}]; femData = state["FiniteElementData"]; femData["PDECoefficientData"]["All"] {{{{0}}, {{{{0}, ...

8

params = {ν1 -> 1, ω1 -> 10, F -> 4}; system = {D[ x1[t], {t, 2}] == -ν1 D[x1[t], t] - ω1^2 x1[t] + F Cos[ω t], x1[0] == 1, x1'[0] == 0}; soln = DSolve[system /. params, x1[t], t][[1, 1]]; (* and the steady state is*) lim = ((List @@ (Expand@soln[[2]])) /. x_ /; (Limit[x, t -> Infinity] == 0) :> 0) steadyState = Simplify@Together[Plus @@ lim] ...

7

With eqn[{k_, r_, H0_, P0_}] := {H'[t] == r (1 - H[t]/k) - d H[t] P[t], P'[t] == -s P[t] + e H[t] P[t], H[0] == H0, P[0] == P0} d = 0.01; s = 0.3; e = 0.02; I would define one simulation as sim := Module[ {k = RandomVariate[NormalDistribution[150, 20]], r = RandomVariate[NormalDistribution[0.4, 0.003]], H0 = RandomVariate[UniformDistribution[{50, ...

7

The curvature and torsion are rates of turning of the Frenet-Serret frame and can be used to integrate the frame using an Euler-type method. The unit tangent vector of the frame is the velocity and can be used to integrate the position. Set up some initial data: the initial Frenet-Serret frame (randomly chosen below), initial point s0, the change in ...

6

Yes! fun = NDSolve[{u ρ'[u] + ρ[u - 1] == 0, ρ[u /; u < 1] == 1}, ρ, {u, 0, 20}, WorkingPrecision -> 40, PrecisionGoal -> ∞, AccuracyGoal -> 40, Method -> "StiffnessSwitching"][[1, 1, 2]]; fun[15] (* 7.5899080042980595046528227779709126741*10^-20 *) DickmanRho[15] (* 7.5899080042980595047*10^-20 *) It works not only in V10 (in ...

5

Use the method option Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> "Projection"}} This forces the equations to be incorporated as constraints. See tutorial/NDSolveDAE#128085219. Depending on the version, you might need to us Rationalize to make the coefficients exact to avoid 1/0 errors. (In general, I avoid machine ...

5

You can use Method -> "FiniteElement" such as: h = 10.6; F = 0.001; d = 1.0; L = 100*d; phi[x_] := Piecewise[{{0.5*(1 - Tanh[x]), x < L*0.5/d}, {0.5*(1 + Tanh[x - L/d]), x > L*0.5/d}}] s = NDSolveValue[{u''[x] == (h)*phi[x]*phi[x]*u[x] - F*d*d*(1 - phi[x]), u[-2.5] == 0.0, u[L + 2.5] == 0.0}, u, {x, -2.5, L + 2.5}, Method -> ...

5

Building on @SjoerdC.deVries answer you can use: sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, DirichletCondition[u[x, y] == 100., x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= 60 || 40 <= x <= 60 && y == 40 || 40 <= x <= 60 && y == 60], u[x, 0] == u[x, 100] == ...

5

Actually a very simple modification would make your approach work: Plot[Evaluate[u[4., x] /. %], {x, 0, 1}] the problem is that for pattern matching (which is what ReplaceAll (/.) does in the background) 4 (with head Integer) and 4. (with head Real) are two different things. You even might meet cases where using 4. would not be enough to make this work, ...

4

From the comments, we can set up the OP's DEs as follows and show they can be solved exactly. First the system is the direct product of two independent systems, so let's separated them. γ = 6; g = -98/10; yIVP = {y''[t] + γ*(y'[t])^2 == g, y[0] == 0, y'[0] == 15/10}; qIVP = {q''[t] == -γ*(q'[t])^2, q[0] == 0, q'[0] == 7}; nysol0 = NDSolve[yIVP, {y}, {t, ...

4

The problem is that NDSolve is not HoldAll or HoldFirst. Therefore the differential equation is evaluated symbolically before it is passed to NDSolve. Thus the differential equation that NDSolve sees is {R*Q'[t] + Q[t]/C2 == newfSign[t], Q[0] == 0} (* {4.54545*10^6 Q[t] + 2000 Derivative[1][Q][t] == 5, Q[0] == 0} *) The reason that one does not see the ...

4

For me it looks like the OddQ function behaves unexpected If it is replaced by OddQ2[n_] := If[Mod[n, 2] == 1, True, False] then s = NDSolve[{R*Q'[t] + Q[t]/C2 == newfSign[t], Q[0] == 0}, Q, {t, 0, 4 \[Tau]}] Plot[Q[t]/C2 /. s, {t, 0, 4 \[Tau]}, PlotRange -> Full] outputs As I stated in my comment, for functions with fast oscillations there are ...

4

NSolve and Reduce can solve this equation by restricting search. As can be seen by inspection: x=0 is a solution and there are two other solutions symmetric about the origin (a plot reveals and noting if r is a root then so is -r: $c \sin(-r)- (-r)= -(c\sin(r)-r)=0$. Quiet@Reduce[y[x] == 0 && Abs@x < 1, x] yields: x == -0.786647 || x == 0 || ...

3

One can do it semi-automatically. Let us introduce a normalized variable $$\xi = \frac{x}{s(t)}, \quad \xi \in [0,1]$$ and make a simple finite difference method over $\xi$. The differential equation in new variables is ClearAll[u, s, x, t, ξ] D[u[x/s[t], t], t] == D[u[x/s[t], t], x, x] /. x -> ξ s[t] If we divide the interval $[0,1]$ by $n$ ...

3

The linear system is easily solved generally by first calulating the MatrixExp. Then we can extract the real and imaginary parts. Here we go The complex matrix A = {{0, 1}, {-2, -I}}; The matrix exp At = MatrixExp[t A] (* {{1/3 (2 Cos[t] + Cos[2 t] + I (2 Sin[t] - Sin[2 t])), 1/3 (I (-Cos[t] + Cos[2 t]) + Sin[t] + Sin[2 t])}, {1/ 3 (2 I (Cos[t] ...

3

It seems to me that using MakeBoxes in this case is overkill. How about this simpler definition? supressVariable[f_Symbol] := Format[f[t, x], TraditionalForm] := Interpretation[f, f[t, x]] SetAttributes[supressVariable, Listable] supressVariable[{v, ρ, p, f}]; This doesn't encounter the issue you faced, because the symbol f is passed directly to ...

3

InterpretationBox holds its arguments (it has HoldAllComplete). You must evaluate ToBoxes[f] outside of this head, easily accomplished with Function as follows: supressVariable[f_Symbol] := f /: MakeBoxes[f[t, x], TraditionalForm] := InterpretationBox[#, f[t, x]] & @ ToBoxes[f]

3

An extended comment that is too long for the comment section. Amplifying on the comment by Stephen Luttrell: eqn = y'[x] + 2 E^x y[x] - y[x]^2 == E^(2 x) + E^x; The generic solution solnG = DSolve[eqn, y[x], x][[1, 1]] y[x] -> E^x + 1/(-x + C[1]) Verifying that the generic solution satisfies the equation eqn /. NestList[D[#, x] &, solnG, 1] ...

3

As it is (unspecified k[t], A) NDSolve will not work. However the equations can be handled analytically. After a simple manipulation you can decouple them and get : rawx[t_] = x[t] /. DSolve[{k[t] x'''[t] - k'[t] x''[t] == -k[t]^3 x'[t]}, x[t], t] rawy[t_] = y[t] /. First@DSolve[{D[#, {t, 2}]/k[t] == y'[t]}, y[t], t] & /@ rawx[t] Now you can check ...

3

I realize from the idiomatic way this program was written that Mathematica is not a language you're very familiar with yet, at least not with the efficient and readable ways Mathematica can address your project. I won't try to suggest every possible improvement but focus on the specific problem of getting the functions you want into the Plot. The following ...

3

I could not verify the solution given by M either. Maple solved this and verifies the solution. But the solution is given as implicit. Here is the solution fyi in case it might help see what is the problem: restart; eq:=diff(u(t),t)=sqrt(u(t))+1/(n+1); ic:=u(0)=0; sol:=dsolve({eq,ic},u(t)); DEtools[remove_RootOf](sol); ...

2

To understand what's going on in this problem, note that $$x'' = - (x + y)\\ y'' = -(x+y)$$ Then just define two new variables describing the sum and difference: $$z \equiv x+y\\ r\equiv x-y$$ and calculate their equations of motion by adding and subtracting the first equations: $$z'' = -2 z\\ r'' = 0.$$ This is a decoupled system of equations, and ...

2

You can get a path function directly from the solution of your ODEs. I don't understand your ODE's, so I'm going to work with a much simpler system, which gives the path of particle moving under constant gravity in a vacuum. g = -9.8; numSoln = NDSolve[{ y''[t] == g, q''[t] == 0., y[0] == 0., q[0] == 0, y'[0] == 50., q'[0] == ...

2

The condition x[t] == 0 && p[t] > 0 is not true very often, if at all, for the initial conditions tried in the OP's differential equation. It would not explain why your session would encounter a such a problem. (I haven't gotten such an error in a very long time. I always associated it with a bug. It was often unreproducible. If you can ...

2

There are at least two approaches by which you could obtain a closed form emulation of your answer. Both involve extracting a list of points that are part of the solution. 1. Inexpensive but takes you outside Mathematica. Export the list of points to a CSV file. Obtain (free trial) a program called Eureqa (http://www.nutonian.com/products/eureqa/) that ...

2

The reasons for the change in the behavior of ParallelTable are subtle. The main source of the problem is that in funcB, the argument k_ is not protected with ?NumericQ like this: funcB[t_?NumericQ, k_?NumericQ] := (* a solution *) funcB[t, k] = Exp[NIntegrate[funcA[et, k], {et, tini, t}]] But more on that later. The problem does not appear in the ...

2

There is a problem with boundary conditions. Changing them and fixing a few typos you can get: eq = {-f[x]^2 + f[x]^4 + x^4 Derivative[1][f][x]^2 + (x^3) /(-1 + x) f[x] ((1 - 2 x) f'[x] - (-1 + x) x f''[x]) == 0, f[1/10] == 1, f'[1/10] == -1/10}; eq // Column // TraditionalForm s = NDSolve[eq, f, {x, 1/10, 9/10}]; Plot[Evaluate[f[x] /. ...

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