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


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

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


12

We can use FEM as it is to solve this problem in versions 12-12.1.1 Needs["NDSolve`FEM`"]; k = 2; Pr0 = 1; Ra = 10^4; R = Ra*Pr0; a = 0; d = 1/2; b = 1/3; Da = 1; Ha = 1; Q = 1; reg = Rectangle[{0., 0.}, {1., 1.}]; {UX, VY, PK, TK} = NDSolveValue[ {{Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}} . Inactive[Grad][u[x, y], {x, y}], {x, y}] + ...


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


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


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


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

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

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

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


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


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