34

Okay, this is a bit of an embarassment. Here is a very small modification of the original code. I simply made explicit option settings, made a denominator to Sin explicitly real, that kind of thing. My tests show the same timing as the original, give or take an iota. ie = 200; ez = ConstantArray[0., {ie + 1}]; hy = ConstantArray[0., {ie}]; fdtd1d = ...


18

Since a fully NDSolve-based solution is acceptable for you, let me give you one. You simply need the magic of "Pseudospectral" or a dense enough 4th order spatial discretization: mol[n_, o_:"Pseudospectral"] := {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n, "DifferenceOrder" -> o}} ...


14

The last column seems in error. Here's a workaround for the sample problem, although it does not fix NDSolve`FiniteDifferenceDerivative: dx2 = dx1; dx2[[All, -1]] = -Reverse@dx1[[All, 1]] (* {-0.50000000000000000000, 0.25989153247414500869, -0.29289321881345247560, 0.36161567304292239214, -0.50000000000000000000, 0.80995720221088751026, -1....


12

Introduction A lack of time to write up an answer ironically provided time to reflect on the problem, and some nagging uncertainties about some issues contributed to the delay. The slowness of the OP's code can be seen by a simple analysis, which reveals that some expensive calculations are repeated multiple times for each step in the time integration. A ...


12

As mentioned in the comment above, handling irregular region with FDM is cumbersome and frustrating in my view, actually that's where I stopped my self-learning of FDM and turned to finite element method (FEM), which is more suitable for this task, and finally write this. (Have a look at this post for more information. ) Nevertheless, there seems to be no ...


12

Let's Compile. dif = 1/500; pec = 20; U0 = 0; V0 = 0; F0 = 3; dn0 = 1; kap = 5; n = 63(*31*); n1 = n + 1;(*dt=40./n^2;*)sm = 4 200; r = 20;(*a=dt dif n n;*) (*c=ConstantArray[0.,{n1,n1}];*) (*d=ConstantArray[0,{n1,n1}];*) den = ConstantArray[dn0 (1 + .3 Tanh[-kap Range[-n1/2, n1/2]]), n1]; (*c0=ConstantArray[0,{n1,n1}];*) u0 = ConstantArray[0., {n1, n1}]; ...


11

The periodic driving at one point doesn't seem to be compatible with the boundary conditions expected by NDSolve, so I modified the problem in two ways: first, broaden the point source into a Gaussian, and then incorporate this driving as a source term in the actual differential equation. So we're actually solving the inhomogeneous wave equation here. For ...


11

The slowness is due to the fact that several steps in your code were not compilable because they invoked MainEvaluate. I localized all the variables by adding Module. Then, several variables were mis-recognized as integer when they should be reals. To fix this I added decimal points to some numbers like 1. The externally defined functions at the beginning ...


11

Partial ND It's possible to again extract all this weight data and stuff in the higher dimensional cases using effectively the same procedure as for the 1D: fdmWeightDataND~SetAttributes~HoldAll; fdmWeightDataND[a_, dep_Symbol, indeps : {__Symbol}, hs : {__Symbol}] := Block[{dep}, Block[indeps, Block[hs, Module[{weights, order, denomWeight, denom, ...


10

NDSolve-based solution To solve the equation set with NDSolve, we need to resolve several issues: As mentioned by bbgodfrey in the comment above, Abs can't be differentiated properly in Mathematica. (This design is reasonable: do remember argument of Abs can be a complex number! ) This can be circumvented by rewriting Abs with Sqrt[(*……*)^2] i.e. modifying ...


10

First of all, if you just want to solve the equation, NDSolve is enough: L = 6; tfin = 1; vmax = 12; ρmax = 200; With[{rho = rho[t, x]}, eq = D[rho, t] + vmax (1 - 2 rho/ρmax) D[rho, x] == 0; ic = rho == Piecewise[{{ρmax, x <= L/2}}] /. t -> 0; bc = {rho == ρmax /. x -> 0}]; mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines", "...


9

To illustrate the problem, I will give an example that differs from the one proposed by Riku. But in this case, numerical instability is better seen. The result is similar to erosion. Perhaps geologists will like this. b = {1, HeavisideTheta[x - y]}; L = 4; reg = DiscretizeRegion[Rectangle[{-L, -L}, {L, L}], MaxCellMeasure -> .01]; eq = D[u[t, x, y], t]...


8

From the tutorial Line Search Methods, there is an example similar to this: Newton's method effectively uses a quadratic model and solves the equation $$H {\bar s} = - \nabla f(x,y)\,,$$ where $H$ is the Hessian $H = \nabla^2f(x,y)$, for the step $\bar s=(\Delta x, \Delta y)$. For an objective function that is a quadratic function like x^2 + 10 y^2, this ...


8

Edit: I have reorganized my earlier answer and added a significant amount of new material. Transformed Equations As suggested in the question, the equations are simpler, if A and B are replaced by their reciprocals. eq1 = (μ^2 B[r] + 1) ϕ[r]^2 + A[r] ϕ'[r]^2 + λ ϕ[r]^4/2 + A'[r]/r + (A[r] - 1)/r^2; eq2 = (μ^2 B[r] - 1) ϕ[r]^2 + A[r] ϕ'[r]^2 - λ ϕ[r]^4/2 ...


8

Here is a code I wrote a while back to get you started: Needs["NDSolve`FEM`"] FiniteElementDerivative[order : {__Integer}, mesh_ElementMesh] /; 1 <= Length[order] <= 3 := Block[{dim, nr, vd, sd, mdata, ccoef, pos, dcoef, cdata}, dim = Length[order]; nr = ToNumericalRegion[mesh]; vd = NDSolve`VariableData[{"DependentVariables", "Space"} ->...


8

Here is how to obtain the "$n$-th approximating function of CantorStaircase. ClearAll[f]; f[0] = x \[Function] x; f[n_Integer?Positive] := f[n] = x \[Function] Evaluate[PiecewiseExpand[ Piecewise[{ {0, x <= 0}, {1/2 f[n - 1][3 x], 0 <= x <= 1/3}, {1/2, 1/3 <= x <= 2/3}, {1/2 + 1/2 f[n - 1][3 x - 2], 2/3 &...


7

NDSolve has trouble in handling the last b.c., so let's help it a bit by discretizing the PDE and corresponding i.c. and b.c. to a set of DAE. First, interpreting your equation set to Mathematica code: γ = 1/100; α = 1; κ = 1; R = 10; Z = 1; eps = 10^-1; tend = 1; eq = With[{u = u[r, Sqrt[κ] z, t]}, Laplacian[u, {r, th, z}, "Cylindrical"] == D[u, t] /...


7

Plugging the solutions into the PDE yields for soltraditional (I D[u[t, x], t] + D[u[t, x], {x, 2}] - I Sin[x] u[t, x]) /. u -> soltraditional; Plot3D[Evaluate@ReIm@%, {x, -Pi, Pi}, {t, 0, tend}, PlotRange -> All, ImageSize -> Large, AxesLabel -> {x, t, u}, LabelStyle -> {Bold, Black, 15}] which is not so good, the spiky behavior near t ...


7

I test nonlinear FEM using solutions obtained by other methods. I developed one of these methods for the problem of the nature convection, aerodynamics, and unsteady hydrodynamics, using linear FEM - see Solver for unsteady flow with the use of Mathematica FEM . Let me give an example. After the release of Mathematica version 12, I tested a non-linear FEM on ...


7

OK, let me extend my comments to an answer. Your code doesn't give proper result because you haven't removed the redundant equations properly. First of all, notice pdetoae will discretize equations in the following way: If the equation is defined on the whole domain of definition, difference equations will be generated on every grid points. In your case, ...


7

It's because NDSolve is using different method for these 2 problems. For the first problem, Shooting method, which is for nonlinear boundary value problem of ordinaray differential equation, is used; while for the second problem, FiniteElement method, which supports nonlinear problem as of v12, is used. Different solvers of NDSolve are not equally powerful. ...


7

Introduction This uses implicit finite difference method. Using standard centered difference scheme for both time and space. To make it more general, this solves $u_{tt} = c^2 u_{xx}$ for any initial and boundary conditions and any wave speed $c$. It also shows the Mathematica solution (in blue) to compare against the FDM solution in red (with the dots on ...


7

NDSolve currently can't handle coupled PDE and ODE, so let's discretize the system all by ourselves, but before that, I'd like to point out the system actually has two solutions. By eliminating $\frac{∂f}{∂r}$ from $(3)$ and $(4)$ we obtain: c = 10; d = 1; e = 1; With[{f = f[r, t], g = g[t]}, bc = {D[f, r] == c D[g, t], f - (d - e D[f, r]) D[f, r] == g} /....


7

If you use the Finite Element Method, no flux is the default boundary condition, so there is no need to specify. An alternative to Daniel's answer would be: (* Define parameters *) l = 6; tend = 0.1; parms = {d -> 2, da -> 5.5, h -> 0.5, k -> 0.5, x0 -> 0.2}; (* Create Parametric PDE operators for n and a *) parmnop = D[n[t, x], t] - d D[n[...


7

Something like the following should do. It employ the finite element method. Ω = DiscretizeRegion[Rectangle[{-10, -10}, {10, 10}], MaxCellMeasure -> (1 -> 0.5)]; sol = NDSolveValue[ { D[Ψ[x, y, t], t] == I/2 Laplacian[Ψ[x, y, t], {x, y}] - I ((x^2 + y^2) + (x + y) Sin[t]^2) Ψ[x, y, t], DirichletCondition[Ψ[x, y, t] == 0, True], Ψ[x, y,...


7

We can use the original code after small correction as follows Clear[gab, \[Rho]0, \[HBar], \[Omega]0, mb, \[Tau], k, r0, t, s, drc] drc = 0; \[HBar] = 1; \[Omega]1 = 10; \[Omega]0 = 10; gab = 1; mb = 1; ma = 1; \[Rho]0 = 1; drc = 0; r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)]; f[\[Tau]_] := (Sqrt[ 2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/ 2) r0^2 \[...


6

This is clearly a bug in 10.0 up to 10.2. According to reply to my report [CASE:3484187] The issue has been resolved in version 10.3. Please upgrade. I hope this will be helpful to someone.


6

You could try using DifferenceDelta to check your answers for these examples. In[1]:= DifferenceDelta[f[x], {x, 0}] Out[1]= f[x] In[2]:= DifferenceDelta[f[x], x] Out[2]= -f[x] + f[1 + x] In[3]:= DifferenceDelta[f[x], {x, 2}] Out[3]= f[x] - 2 f[1 + x] + f[2 + x] In[4]:= DifferenceDelta[f[x, y], {x, 0}, {y, 1}] Out[4]= -f[x, y] + f[x, 1 + y] In[5]:= ...


6

You can use FindRoot to find zeros of the derivative. FindRoot as an option DampingFactor that should serve your purpose: ListLinePlot[ Last[Reap[ FindRoot[{D[x^2 + y^4, x] == 0, D[x^2 + y^4, y] == 0}, {{x, 10}, {y, 10}}, DampingFactor -> .1, EvaluationMonitor :> Sow[{x, y}]]]], PlotRange -> All, PlotMarkers -> {Automatic, 10}] ...


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