New answers tagged

0

If you solve for the second derivatives, you won't have to use "EquationSimplification" -> "Residual" and things will work ok. Solving for the second derivatives be faster if you start with exact coefficients. Also, if you clear l, solve for the derivatives, and then substitute a value for l, Solve won't choke on the algebra. The long time it takes is ...


2

I have experienced serious memory leaks with NDSolve, NIntegrate, and FindRoot in every version of Mathematica I have owned (still on V9), perhaps over 20 years. They typically show up for me when the routine is called at a deep level in a complex program. The only cure in many cases is to remove the offending routine and write your own. For NIntegrate, ...


3

This turns out to be a bug since v9. v8.0.4 gives the correct result: $Version eqn = D[x[t], t] == a*x[t] + c*t^n; bc = x[0] == x0; sol = Simplify[DSolve[{eqn, bc}, x[t], t], {a > 0, n >= 1, c > 0}] (* "8.0 for Microsoft Windows (64-bit) (October 24, 2011)" *) (* {{x[t] -> a^(-1 - n) E^(a t) (a^(1 + n) x0 + c Gamma[1 + n] - c Gamma[1 + n, a t])}}...


1

An alternative approach is to solve for t[x] instead of x[t]. t[x] /. First@DSolve[t'[x] == 1/(5 x^4 + 3 x^(-4/3)), t[x], x] /. C[1] -> 0; ParametricPlot[{Chop@%, x}, {x, 0, 3}, AxesLabel -> {t + C[1], x}, AspectRatio -> GoldenRatio] Because the ODE determines t[x] only up to an arbitrary constant, the curve above can be shifted by an ...


2

Making my comment an answer: Try to extract the mesh from the eigenfunction and use that for integration. Something like: NIntegrate[ circlepluck[x, y]*circlefuncs[[i]][x, y], {x, y} \[Element] circlefuncs[[1]]["ElementMesh"]] This would switch off the adaptive mesh refinement.


6

Since a fully NDSolve-based solution is acceptable for you, let me give you one. You simply need the magic of "Pseudospectral": mol[n_, o_:"Pseudospectral"] := {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n, "DifferenceOrder" -> o}} ic[θ_] = Piecewise[{{1.001, -2 Pi/3 < θ < -...


0

There is no library that I am aware of that does the trick. So there is no easy solution on it. Anyway, there is a publication about CUDA and parabolic PDEs. Do a search on "CUDA Jahn Uni Bayreuth" in the internet. There is a pdf (diploma Dissertation). May be that helps to find a solution for your problem; adding the problem to learn German ;).


7

Alternatively, you can use Needs["NDSolve`FEM`"] \[CapitalOmega] = ToElementMesh[ RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]], MaxCellMeasure -> {"Area" -> 1}] This will generate a second order mesh and you'll need far fewer elements to get an accurate solution then with a first order mesh (as ...


3

Using DiscretizeRegion[] with RegionDifference[] prior to NDSolveValue was the key. Ω = DiscretizeRegion[ RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]], Method -> "Continuation", AccuracyGoal -> 7, PrecisionGoal -> 7, MaxCellMeasure -> {"Area" -> 0.1}]; sol = NDSolveValue[{D[u[x, y], x, x] + ...


2

ClearAll["Global`*"] Remove["Global`*"] {ysol, xsol} = ParametricNDSolve[{x'[t] == x[t] + g*x[t]*y[t], y'[t] == 1 - 2*x[t]^2 - g*y[t]^2, x[1] == 1, y[1] == 1}, {y, x}, {t, 0, 10}, {g}]; ParametricPlot[Evaluate@Table[{y[g][t] /. ysol, x[g][t] /. xsol}, {g, 0, 1, 1/4}], {t, 0, 2}, PlotRange -> {{-3, 3}, {-3, 3}}, PlotLegends -> {"g=1/4", "g=2/4", "g=...


2

Maybe like this: sol = DSolve[x'[t] == 5 x[t]^4 + 3 x[t]^(-4/3), x[t], t] eq = sol[[1, 1, 2, 0, 1]][x] == (sol[[1, 1, 2, 1]] /. C[1] -> 2) ContourPlot[sol[[1, 1, 2, 0, 1]][x] == (sol[[1, 1, 2, 1]] /. C[1] -> 2), {t, -3, 3}, {x, -3, 3}, Axes -> True, Frame -> False, AxesLabel -> {t, x[t]}] EDITED: If You exectue this code: Internal`...


2

Your initial and boundary conditions are inconsistent. Check NDSolve::ibcinc how to avoid that. Also I think you need additional boundary condition. Anyway I'm not sure your constant functions are good for boundary/initial conditions. The r domain includes zero and it gives error in your $1/r$ term. Using the recipe from above link and adding another ...


8

Setup The region to solve the PDEs Clear["Global`*"] xMax = 453.595 - Sqrt[450.05^2 - 3^2]; regA = Rectangle[{0, 0}, {xMax, 6}]; regB = ImplicitRegion[{453.595 - Sqrt[450.05^2 - (y - 3)^2] - x < 0}, {{x, 0, xMax}, {y, 0, 6}}]; reg = RegionDifference[regA, regB]; RegionPlot[reg, AspectRatio -> 1/10] Solve the PDEs using NDSolve solV = ...


4

One can use the new function DistanceMatrix[] for the purpose; this avoids repeated computations (since the underlying matrix is symmetric). With[{stepSize = 2, end = TMax}, MatrixPlot[UnitStep[0.01 - DistanceMatrix[uSolpbc[Range[0, end, stepSize], 0]]]]] and your plot is produced very quickly, without the need to invoke parallelization.


2

The plot in the question must have been obtained with stepSize = 15, not stepSize = 2. Using the latter value gives a smooth plot, The computation takes about 78 sec on my PC. To address the specific issue in the question, the run time can be reduced by two orders of magnitude using Block[{stepSize = 2, end = TMax, tt, rd}, tSolpbc = Table[uSolpbc[...


2

No, that's not possible. If you could show your code there might be other ways to speed it up.


1

Your problem is twofold. The first part is very subtle and sneaky. Your "c" characters are actually typed in as the Unicode Character 'cyrillic small letter es' (U+0441) in the NDSolve expression, as you can see from the following: Inactive[NDSolve][{ I*с00'[t] == с10[t] + с01[t] + с00[t] + с00[t], I*с01'[t] == с11[t] + с00[t] + с01[t] - с01[t], ...


1

Thank for the help with @xzczd ,I want to post a possible answer. The code is kp = 2; vm = 311; w = 100*Pi; k = 1; ki = 100; nsol = NDSolveValue[{y''[ x] + (kp*vm*Cos[w*x - y[x]] - 2/3 kp*k*vm*Cos[w*x] Cos[y[x]]) y'[x] - ki*vm*Sin[w*x - y[x]] - 2/3 ki*k*vm*Cos[w*x] Sin[y[x]] - kp*vm*w*Cos[w*x - y[x]] == 0, y[0] == 0, y'[0]...


0

After digging around for a bit I found the tutorials in the standard documentation: FEMDocumentation/tutorial/SolvingPDEwithFEM and FEMDocumentation/tutorial/FiniteElementProgramming. Both have sections on methods to trade memory/time/accuracy, but the first tutorial is way easier for people new to Mathematica! The most promissing proposals seem to be: ...


2

Let's read the error message: NDSolve`Iterate::ndsz: At t == -1.*10^10, step size is effectively zero; singularity or stiff system suspected. The "step size is effectively zero" means that in floating-point arithmetic t + dt is equal to t for the computed time step dt. If this is the problem, increasing working precision might help. de = {F'[t] == -(...


3

If you accept that the general solution can be constructed as an (infinite) Fourier series, $$\sum_{m=-\infty}^{\infty}\phi_m(r) \exp(i m \theta),$$ then you can obtain the expression for $\phi_m(r)$ as follows: eqn = 0 == Simplify@ Laplacian[ Laplacian[ ϕ[r] Exp[I m θ], {r, θ}, "Polar"], {r, θ}, "Polar"]; DSolve[eqn, ϕ, r] (* ==> {{ϕ -> ...


0

for example for function y'=f(x,y)=x-y , I think the easiest way is... f[x_, y_] := x - y yn[x_,n_,x0_,y0_]:=If[n>0, y0+Integrate[f[s,yn[s,n-1,x0,y0]],{s,x0,x}], y0] example (y(0)=1 and 3 iterations): yn[x, 3, 0, 1] 1 - x + x^2 - x^3/3 + x^4/24


5

Maybe commandeer FindRoot with the Villegas-Gayley trick: Updated, with the order of the steps taken by FindRoot saved in icsteps. The results of FindRoot, as saved by NDSolve and shown below as DownValues[], have been sorted by Mathematica and are not in the order in which there were called. This update stores the order in icsteps. Clear[x, t]; Internal`...


4

I think this gives what you want: We construct the integral inside the definition of H by adding another ODE to the NDSolve system, which I called logH. This in fact calculates the integral from ic, not from 0. So to define H we need to subtract logH[0] from logH[t] before exponentiating. This should be a much more accurate (and faster) way of computing ...


4

As you have noted placing the parameter on the left hand side of a boundary condition: y[b] == 1 does not work. You are forced to keep the parameter on the right hand side. y[0] == a What can be done is to keep the parameter on the RHS and then use numerical methods to determine the value of b for a particular a parameter (in other words, seek a ...


1

The problem with your system of equation is that the "zero" far field solution is unstable, hence extremely sensitive to the initial conditions. Posing the problem as an initial value problem, with the "known" conditions from the successful solution: max = 50; Pr = .72; a = 0.7172594734816521` b = -0.4344414944896132` pohl = NDSolve[{f'''[\[Eta]] + 3 f[\[...


2

Here is a function for multiplying InterpolatingFunctions generated from a single NDSolve (so that the coordinate grids are the same, as well as one-dimensional). This yields a single InterpolatingFunction that interpolates the product. It carries over derivative information, too. (I've done this before on the site, I think.) Anyway, Integrate on an ...


3

You have many superfluous sets of {} that generate unexpected output in your code. In particular, the Interpolation functions generated by NDSolve were not returning a scalar value, but instead a unidimensional vector, i.e. a list containing a single value instead. That was probably an unintended consequence of the extra sets of braces in the definitions of ...


3

Clear["Global`*"] yy = 10^-4; rr = 0.999; xx = 10^-15; zz = 10^-4; mm = 10^-4; yy + rr + xx^2 + zz - mm^2 - zz^2/24 ic = -17.5 s = NDSolve[{D[y[t], t] == (3 y[t])/5 - (12 m[t]^2 y[t])/5 + (2 r[t] y[t])/ 5 - (6 x[t]^2 y[t])/5 + (3 y[t]^2)/5 + (7 y[t] z[t])/ 5 - (y[t] z[t]^2)/10, D[r[t], t] == -((2 r[t])/5) - (12 m[t]^2 r[t])/5 + (2 r[...


7

The problem is with the default starting initial conditions used by the shooting method in NDSolve. The shooting method is where FindRoot is being used internally, so the OP's error message is a strong hint that this is the problem. Getting convergence in a nonlinear system can depend greatly on the starting conditions. Having luckily solved the system ...


2

The singularity at r == 0 is of the type $O(1/\sqrt{r-1})$ and is easily removed with a substitution $r - 1 = u^2$. As for $r = \infty$, f[r] has a bunch of factors of the form $$ (A + B r^n)/(C r^n)$$ plus the final two factors $(r-1)r$. This means that the function is asymptotic to $a r^2$ at $r = \infty$ for some constant $a$. This makes the integrand ...


3

The e in the non-linear fit should be Exp[] or it needs to be a variable x = Clip[r[t], {0, 1}]; y = Clip[rpp[t], {0, 1}]; tdata = NDSolve[{r'[t] == -L1*rpp[t] + c, rpp'[t] == L2*rpp[t] (1 - (rpp[t]/(r[t] + rpp[t]))) + L3*r[t], r[0] == 1, rpp[0] == 0} /. {L1 -> 0.5, c -> 0.001, L2 -> 0.91, L3 -> 0.05}, {r[t], rpp[t]}, {t, 0, 10}...


0

As mentioned in the comments above, your understanding for function is wrong. You'd better put more effort in learning the core language of Mathematica, or even consider revisiting the definition of function in math. The following is a not-the-simplest but correct way to code your equation set in Mathematica, do read it carefully: Clear["`*"] vei = 1/1000; ...


1

You cannot use square brackets for parentheses and must include initial conditions. eqns = { x'[t] == -z[t] (x[t] + 1.5 Log[z[t]]), y'[t] == z[t] - 1, z'[t] == -(.64) (y[t] + (x[t] + Log[z[t]])/.8), x[0] == y[0] == z[0] == 1} // Rationalize // Simplify; sol = NDSolveValue[eqns, {x[t], y[t], z[t]}, {t, 0, 1000}]; Manipulate[ ...


2

This works Clear[x, y, z] sol = NDSolve[Join[ Thread[{x'[t], y'[t], z'[t]} == {-z[t] x[t] + 1.5 Log[z[t]], z[t] - 1, -(0.64) y[t] + (x[t] + Log[z[t]])/0.8}], Thread[{x[0], y[0], z[0]} == {1, 1, 1}]], {x, y, z}, {t, 0, 100}]; ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 100}]


5

As @yarchik pointed out, this is strictly speaking not a two-electron problem, but a problem of "distinguishable particles," since the antisymmetry postulate (and spin) is ignored. Moreover, the 1D Coulomb potential in the initial equations can't be just $1/x$ -- it has to have an absolute sign because otherwise it would be attractive on one side and ...


0

I have found yet another way to solve this. Based on J.M.'s answer and another post on defining the derivative of Abs for Real valued functions (will look it up later to add here). The following code works: Derivative[1][Abs][x_] := Re[ Conjugate[x] D[x, s] ]/Abs[x]/D[x, s]; Plot[D[Abs[k[s]], s] /. k -> kf /. s -> ss, {ss, 0, 1}] It is worth ...


1

When you see the warning NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help. This may indicate the NDSolve did try to solve this as a purely spatial problem. You can force it so use the MethodOfLines via the option Method -> "MethodOfLines" for time integration. This then solves in ...


2

Alternatively, you could do a little complex number algebra for the purpose: kf = NDSolveValue[{I k'[s] == 20 (s k[s] - (1 - s) (1 - k[s]) (1 + k[s])), k[0] == 1/2}, k, {s, 0, 1}]; Plot[{10 Abs[kf[t]], Re[kf[t] Conjugate[kf'[t]]]/Abs[kf[t]]}, {t, 0, 1}, PlotLegends -> {"scaled function", "numerical derivative"}]


1

I may have found a nice way. However, I think the output is a bit ugly. Here's the working code : NullCurve1[n_] := NDSolve[{ u'[r] == 2/(1 - 1/r), u[0] == 2 n }, {u}, {r, 0, 1}, Method -> Automatic, MaxSteps -> 100000 ] NullCurve2[n_] := NDSolve[{ u'[r] == 2/(1 - 1/r), u[5] == 2 n }, {u}, {r, 1, 10}, Method -&...


2

As discussed in meta, I will expand on @J.M.'s comment to provide this question with an answer. Since I dont have the code for your NDSolve I came up with (based on Mathematica documentation): sol = NDSolve[{D[h[t, x], t] == D[h[t, x], x, x], h[0, x] == 3/2, h[t, 0] == 3/2, h[t, 5] == 3/2}, h, {t, 0, 10}, {x, 0, 5}]; This done, you should obtain the ...


2

I will use Eq. 22.7 of your reference $ v = 2 r + 4 m \ln(|r - 2 m|) + B$ The obliques are, I guess is given by $v=-u$. For scaling I am considering $u=2mr$. m = 1/2; Show[ContourPlot[Evaluate@Table[ v == 2 r + 4 m Log[Abs[r - 2 m]] + n, {n, -5, 5}], {r,0,4 m}, {v,-3m,3m}, ContourStyle -> Blue, AspectRatio -> 3/2], ContourPlot[Evaluate@Table[v == -...


1

e = .16; g = .4; w = .97; You must add a initial condition, assume: L[0]=1 sol = NDSolveValue[{y1'[t] == y2[t], y2'[t] == -((2/L[t]) L'[t] + g L[t]) y2[t] - (1/L[t]) Sin[y1[t]], L'[t] == 7 e w Sin[w t + 9 Pi/8]^6* Cos[w t + 9 Pi/8], y1[0] == -1, y2[0] == 1, L[0] == 1}, {y1, y2, L}, {t, 0, 30}]; ParametricPlot[{sol[[1]][t], sol[[2]][t]}, {t, 0, 30}...


4

Preamble After solving the differential equations it appears that there are five equations with five unknowns (four generated from the two differential equations and the parameter k. However there are only four pieces of information in the five equations so there still is one arbitrary constant. There are many routes that lead to this conclusion, ...


3

You could try using numerical differentiation from the Numerical Calculus package: Clear[r, fun] Needs["NumericalCalculus`"] r = NDSolve[ {I k'[s] == 20 (s k[s] - (1 - s) (1 - k[s]) (1 + k[s])), k[0] == 1/2}, k, {s, 0, 1} ]; fun[s_] := Abs[(k /. r[[1]])[s]]; Plot[ {10 fun[t], ND[fun[s], s, t, Terms -> 20]}, {t, 0, 1}, PlotLegends -&...


2

I have revised the code to use StreamPlot instead of VectorPlot and NDSolve. The latter command created issues when the domain for the solution extended to negative times. I don't know why, but NDSolve produced integral curves that did not follow the direction field produced by VectorPLot as provided by MarcoB's nice reformulation of my original code. I will ...


3

Your code cannot work because you defined the potential and electric field with incorrect syntax. Here is a corrected version: Clear[x, y, z]; PhiLineSegment[x_, y_, z_] = (segCharge/(4 Pi e0 segLength)) Log[(segLength/2 + z + Sqrt[(x + segOffset)^2 + y^2 + (z + segLength/2)^2])/(-(segLength/2) + z + Sqrt[(x + segOffset)^2 + ...


5

Her is a start: Clear[ψD]; ψD[j_, k_][x_, y_] := Sin[j Pi x] Sin[k Pi y] Clear[ψN]; ψN[j_, k_][x_, y_] := Cos[j Pi x] Cos[k Pi y] ContourPlot[ψD[1, 2][x, y], {x, y} ∈ Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}], PlotPoints -> 100, AspectRatio -> Automatic, PlotRange -> All] ContourPlot[ψN[2, 2][x, y], {x, y} ∈ Polygon[{{0, 0}, {1, 0}, {1,...


5

Here is a workaround suggested by the response I received from WRI: {sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., WhenEvent[Norm[{x'[t]}] < 0.0001, {x[t] -> 1.}]}, {x}, {t, 0, 200}, Method -> {"EquationSimplification" -> "Residual"} ] Plot[x[t] /. sol, {t, 0, 200}] Warning: This option works by converting the system to a DAE, ...


3

As I mentioned in comments, your main problem was the placement of the definition of z; that definition does not need to be re-evaluated every time there is a modification within the body of Manipulate since it's already delayed. The correct placement for such a definition would be in the Initialization code for the Manipulate. Similarly, the definitions of ...



Top 50 recent answers are included