# FEM derivative matrix construction

Assuming $$A$$ is a $$3 \times 3$$ matrix (and a function of $$x$$, and $$z$$) and $$K_i, \beta_i$$, and $$\alpha(T)$$ are known parameters, I need to solve the following equation, (with an implied sum over $$\mu, i$$, and $$j$$),

$$K_1 \partial^2_j A_{\mu i} + K_{23} \partial_i (\partial_j A_{\mu j}) = 2\beta_1 Tr(A A^T) A^*_{\mu i} + 2\beta_2 Tr(A A^\dagger) A_{\mu i} + 2 \beta_3 [A A^T A^*]_{\mu i} + 2 \beta_4 [A A^\dagger A]_{\mu i} + 2 \beta_5 [A^* A^T A]_{\mu i} + \alpha(T) Tr(A A^\dagger)$$

This system of equations can be simplified with a known form of $$A$$. For example,

$$A = \begin{pmatrix} A_{xx} & 0 & 0 \\ 0 & A_{yy} & 0 \\ 0 & 0 & A_{zz} \end{pmatrix}$$

the equations reduce to (when normalized),

$$a_{||}'' = -\frac{1}{5} a_{||} ( 2 a_{||}^2 + a_\perp^2 ) + \frac{2}{5} a_{||} ( 2 |a_{||}|^2 + |a_\perp|^2 ) + \frac{2}{5} |a_{||}|^2 a_{||} - a_{||}$$

and

$$3 a_\perp'' = -\frac{1}{5} a_\perp^* ( 2 a_{||}^2 + a_\perp^2 ) + \frac{2}{5} a_\perp ( 2 |a_{||}|^2 + |a_\perp|^2 ) + \frac{2}{5} |a_\perp|^2 a_\perp - a_\perp$$

These appear rather simple, and I figured that Mathematica could solve them with NDSolve, but it hasn't worked. With more complicated forms of $$A$$ and more complicated regions, I'm looking to use Mathematica's FEM solver, but I just don't understand how to build the derivative matrix. I've looked at Wolfram's PDE solving guide and their FEM guide, but I only see coefficients for gradients, divergences, etc., while I will also need mixed derivatives ($$\partial_i (\partial_j A_{\mu j})$$). How do I build the mixed derivative matrix in InitializePDECoefficients? I have some working code in C++ that using the finite difference method, and I build the matrices there, but shouldn't there be a better way than inserting every element? Here it seems that the matrices were large, but Mathematica documentation here, they use a $$2 \times 2$$ identity matrix, but for the coefficient? I'm a little confused there.

Or, how can I get NDSolve or NDSolveValue to work? The boundary conditions seem to be causing problems in these functions. The BC's for my equations are, (for the first BC, I just assume that zMax (in code below) is close enough to infinity.)

$$\partial_z A_{\alpha j}|_{z=0} = \frac{1}{b}A_{\alpha j}|_{z=0}, \text{ for } b \in [0,\infty)$$

and

$$\lim_{z \to \infty} A_{\alpha j} = 1.$$

My code so far (using NDSolve first),

b = 1;
{xMin, xMax} = {-5, 5};
{zMin, zMax} = {0, 20};
sol = NDSolve[
{D[p[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[p[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) p[x, z] + 2/5 Abs[p[x, z]]^2 p[x, z] - p[x, z],
3 D[s[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[s[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) s[x, z] + 2/5 Abs[s[x, z]]^2 s[x, z] - s[x, z],
(D[s[x, z], z] /. {z -> zMin}) == 1/b s[x, zMin],
s[x, zMax] == 1,
p[x, zMax] == 1,
(D[p[x, z], z] /. {z -> zMin}) == 0 },
{p, s}, {x, xMin, xMax}, {z, zMin, zMax}
]


This code gives the same error twice (I thought I had other errors previously, but I can't seem to recreate them).

NDSolve: The expression s^(0,1)[x, 0] == s given as a spatial boundary condition for the possibly automatically chosen finite element method should not have explicit derivatives of the dependent variables. NeumannValue should be used to specify spatial derivatives at the boundary.

And my other attempt,

Needs["NDSolveFEM"]
b = 1;
{xMin, xMax} = {-5, 5};
{zMin, zMax} = {0, 20};
nRegion = ToNumericalRegion[Rectangle[{xMin, zMin}, {xMax, zMax}]];

Γpi = NeumannValue[0, z == 0];
Γpf = DirichletCondition[p[x, z] == 1, z == zMax];
Γsi = NeumannValue[s[x, z]/b, z == 0];
Γsf = DirichletCondition[s[x, z] == 1, z == zMax];

eqP = -D[p[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[p[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) p[x, z] + 2/5 Abs[p[x, z]]^2 p[x, z] - p[x, z];

eqS = -3 D[s[x, z], {z, 2}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[s[x, z]] + 2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) s[x, z] + 2/5 Abs[s[x, z]]^2 s[x, z] - s[x, z];

sol = NDSolveValue[{{eqP == Γpi, Γpf}, {eqS == Γsi, Γsf}}, {p, s}, {x, z} ∈ nRegion]


gives these errors, twice:

DiscretizePDE: The FEMStiffnessElements operator failed.

FindRoot: The minimal damping factor of 1/10000 has been reached.

FindRoot: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

NDSolveValue: PDESolve could not find a solution.

• Welcome to MMA Stack Exchange. Please write up your question in this post rather than as an extra link. This will make it much more likely that someone can help you. Dec 31 '21 at 5:01
• @Dunlop Seems that OP is stumbled by the automatic checking. Now the question should be OK. Dec 31 '21 at 5:57
• So, the final system is just an ODE system that doesn't involve derivative of $x$ at all? Dec 31 '21 at 6:07
• @xzczd, thank you for the help editing! I am learning Mathematica's FEM solver, testing it on one of the simplest sets of DEs that come from these equations--the one I describe here does only end up with $z$ dependence, but the next one I want to try will have both $x$ and $z$. Dec 31 '21 at 17:00

I do not fully understand what you want to do but try this:

Needs["NDSolveFEM"]
b = 1;
{xMin, xMax} = {-5, 5};
{zMin, zMax} = {0, 20};
nRegion = ToNumericalRegion[Rectangle[{xMin, zMin}, {xMax, zMax}]];

\[CapitalGamma]pi = NeumannValue[0, z == 0];
\[CapitalGamma]pf = DirichletCondition[p[x, z] == 1, z == zMax];
\[CapitalGamma]si = NeumannValue[s[x, z]/b, z == 0];
\[CapitalGamma]sf = DirichletCondition[s[x, z] == 1, z == zMax];

eqP = -Inactive[Div][{{0, 0}, {0, 1}} .
z}] == -(1/5) (2 p[x, z]^2 + s[x, z]) Conjugate[p[x, z]] +
2/5 (2 Abs[p[x, z]]^2 + Abs[s[x, z]]^2) p[x, z] +
2/5 Abs[p[x, z]]^2 p[x, z] - p[x, z];

eqS = -3 Inactive[Div][{{0, 0}, {0, 1}} .

• Hmmm... it does give a solution! But it's not the one it should be--maybe my BC's don't define a unique solution? Do you know, if I don't give any BC's along the sides (x = 5, x = -5), that it will assume zero? Jan 4 at 0:59
• You could try specifying an InitialSeeding to try to get a different solution. If you do not specify any bcs on some part of the boundary it will implicitly assume a Neumann 0 bc there. Jan 4 at 7:28