4
$\begingroup$

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["NDSolve`FEM`"]
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.

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Dunlop
    Dec 31, 2021 at 5:01
  • 1
    $\begingroup$ @Dunlop Seems that OP is stumbled by the automatic checking. Now the question should be OK. $\endgroup$
    – xzczd
    Dec 31, 2021 at 5:57
  • 1
    $\begingroup$ So, the final system is just an ODE system that doesn't involve derivative of $x$ at all? $\endgroup$
    – xzczd
    Dec 31, 2021 at 6:07
  • $\begingroup$ @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$. $\endgroup$
    – Izek H
    Dec 31, 2021 at 17:00

1 Answer 1

4
$\begingroup$

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

Needs["NDSolve`FEM`"]
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}} . 
      Inactive[Grad][p[x, z], {x, z}], {x, 
      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}} . 
      Inactive[Grad][s[x, z], {x, z}], {x, 
      z}] == -(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 == \[CapitalGamma]pi, \[CapitalGamma]pf}, \
{eqS == \[CapitalGamma]si, \[CapitalGamma]sf}}, {p, 
   s}, {x, z} \[Element] nRegion]

This returns two interpolating functions, not sure if correct.

$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Izek H
    Jan 4, 2022 at 0:59
  • $\begingroup$ 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. $\endgroup$
    – user21
    Jan 4, 2022 at 7:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.