32

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 = ...


26

These days I've picked up some knowledge about finite difference method and now I'm able to fix OP's code :D. Well, before beginning, it should be mentioned that, it's better to avoid For and Subscript in Mathematica, but I'd rather not talk about them in this answer since these have been discussed a lot in this site, what's more, they're not the root of ...


22

The good news is that yes, there is an easy way to put your problem into NDSolve by using the new finite element functionality in v10. The bad news is that it seems the specific problem you're trying to solve is ill-posed. NDSolve can now handle internal boundaries; see e.g. the first figure under "Details" for DirichletCondition. Generating a mesh with ...


17

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}} ...


13

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

Take it in steps: Extract the coefficients and locations into an appropriate data structure. Use that data structure to create the graphics. By examining the FullForm of the original expression, we can cobble a rule to find the key data: the coefficients $c$, $d$, and $-1$ and the offsets to the indexes. First, the expression itself: s = Subscript; exp = ...


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 ...


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

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 ...


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

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, ...


9

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", "...


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

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]...


7

A quick answer now. I will come back to this once I have more time. First we use the common finite difference operators to discretize PDE. Symbolics: Clear[u]; RFSDiscrit[eq_] := Module[{mid}, mid = Distribute@FullSimplify@ExpandAll[( eq /.{ u[x, t] -> u[i, n], D[u[x, t], t] -> (u[i, n] - u[i, n - 1])/\[...


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

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"} ->...


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

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

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, ...


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

Although you haven't exactly asked this, you might like to generate your graphic automatically by applying pattern matching on your difference operator. The basic idea is as below: Clear[i, j, k]; op = Plus @@ MapThread[Subscript[u, i - #1, j - #2, k - #3] &, RotateRight[{0, -1, 1, 0, 0, 0, 0}, #] & /@ {0, 2, 4}] Giving $op = u_{i-1,j,k}+u_{...


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}] ...


5

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]:= ...


5

Just for fun, here is a version using NDSolve: sol = NDSolve[ { x'[t]==-2x[t] Exp[t], y'[t]==-4y[t]^3 Exp[t], x[0]==10, y[0]==10, WhenEvent[x'[t]^2 + y'[t]^2 < 10^-10, end=t; "StopIntegration"] }, {x,y}, {t,0,Infinity} ]; Note that I added a scaling factor (Exp[t]) so that the t range is not too ...


5

As mentioned above, this problem is somewhat similar to, but even more troublesome than this one, because there seems to be no obvious way to discretize and transform the PDE system to an ODE system. To be more specific, it's hard to introduce derivative of $q$ with respect to $t$ without causing other trouble. So, let's give up NDSolve and turn to pure ...


4

This is a bug that has been fixed as of Mathematica 11.0.0.


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