I am trying to solve the following differential equation : $$ \nabla\cdot\left(\begin{pmatrix} c_1 & c_2 \\ -c_2 & c_1 \end{pmatrix}\nabla\phi\right) = 0 $$ on a rectangular domain [1,0]x[0,1] with $c_2$ and $c_1$ constant, and $c_2/c_1\gg 1$.

The boundary conditions are $\phi = \pm 1$ on $x=1$ and $x=0$ respectively, and $n\cdot\left(\begin{pmatrix} c_1 & c_2 \\ -c_2 & c_1 \end{pmatrix}\nabla\phi\right) = 0$ on $y = 0$ and $y=1$.

Here's my confusion :

Since c-matrix is constant, this will evaluate to the Laplace equation, and the Neumann boundary condition will become $\partial_y \phi = 0$, and the information about $c_2$ much larger than $c_1$ is lost, and as expected the contour plot looks like :

sol1 = NDSolve[{D[u[x, y], x, x] + D[u[x, y], y, y] == 
     NeumannValue[0, y == 0] + NeumannValue[0, y == 1], 
    DirichletCondition[u[x, y] == 1, x == 1], 
    DirichletCondition[u[x, y] == -1, x == 0]}, 
   u, {x, y} ∈ 
    ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}]];

ContourPlot[u[x, y] /. sol1, {x, y} ∈ mesh]

enter image description here

But if I set up the equation with a c-matrix then the equipotential lines look strange:

Edit: I forgot to say what omegatimestau is, it’s a function that is proportional to B, in the code snippet below it’s ~ 100, but I get the scratchy lines on the left and right sides of the box even if B = 0.

eqn = Inactive[
     Div][{{1, -omegactimestau[B]}, {omegactimestau[B], 1}} . 
     Inactive[Grad][phi[x, y], {x, y}], {x, y}] == 0;
sol1 = NDSolve[{eqn == 
    NeumannValue[0, y == 0] + NeumannValue[0, y == 1], 
   DirichletCondition[phi[x, y] == 1, x == 1], 
   DirichletCondition[phi[x, y] == -1, x == 0]}, 
  phi, {x, y} ∈ 
   ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}]]

The plot looks like this:

enter image description here

So, I have two questions :

  1. How do I set up the problem so that information in the constant c-matrix is not lost?
  2. What is the strange behavior happening in the second case?

Thank you for reading this.

  • $\begingroup$ What's the definition of omegactimestau[B]? $\endgroup$
    – xzczd
    Feb 26, 2022 at 2:20
  • $\begingroup$ @xzczd omegatimes[B] is a function linear in B it’s a constant in the above operator definition. Sorry I didn’t specify this in the question. $\endgroup$
    – Charlie
    Feb 26, 2022 at 3:02

1 Answer 1


First of all there's a simple mistake in your code, you've included == in the PDE twice! (Once when defining eqn, next inside NDSolve. ) Removing the redundant == results in a reasonable result. Still, it's slightly scratchy at the boundary, so let's make the grid denser:

<< NDSolve`FEM`
omegactimestau[B] = 100;
eqn = Inactive[Div][{{1, -omegactimestau[B]}, {omegactimestau[B], 1}} . 
      Inactive[Grad][phi[x, y], {x, y}], {x, y}] == 0;
reg = ToElementMesh[ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}], 
                    MaxBoundaryCellMeasure -> 10^-3];

sol1 = NDSolveValue[{eqn, phi[1, y] == 1, phi[0, y] == -1}, 
                    phi, {x, y} ∈ reg];   

ContourPlot[sol1[x, y], {x, 0, 1}, {y, 0, 1}, PlotPoints -> 100]

Mathematica graphics

I've omitted the zero NeumannValues, because this is the default setting of FiniteElement method.

  • $\begingroup$ @xzcd Thank you for this answer, I don't understand the pdetoae function very well in the linked answer above, I have two doubts: 1) If I run the above code, it gives me error/warning ? saying 'Part 1 of {}' does not exist and 'Symbol called with 0 arguments'. 2) Do you know how I should go about solving such a differential equation because this is a well defined problem in the physical sense, by that I mean the boundary conditions have physical meaning of a constant voltage on the left and right sides of the box and current is zero on the top and bottom, and the c-matrix is conductivity. $\endgroup$
    – Charlie
    Feb 26, 2022 at 21:54
  • $\begingroup$ @Shubham 1. Have you executed the whole definition of pdetoae? Notice you need to execute the whole code snippet in the link. 2. If this is a well-studied problem, can you add a reference? $\endgroup$
    – xzczd
    Feb 27, 2022 at 2:13
  • $\begingroup$ @xzcd yes, I executed the code snippet in the link of the definition of pdetoae. The problem is well studied because the off-diagonal elements in the c-matrix are the conductivity for when voltage is applied in say, x and current is measured in y which is the Quantum Hall effect. Here's a reference, journals.aps.org/prb/abstract/10.1103/PhysRevB.47.15727. $\endgroup$
    – Charlie
    Feb 27, 2022 at 2:20
  • $\begingroup$ @Shubham Then I guess you don't execute definition of eqn. My eqn is the same as yours so it's omitted in the answer. $\endgroup$
    – xzczd
    Feb 27, 2022 at 2:24
  • $\begingroup$ @xzcd yup, it works. I think there were some conflicts in the code which gave those errors (because I had executed everything before). Anyways, thanks again for pointing out pdetoae, and the answer, as a newbie I am learning a lot. Do you have any advice for solving the pde in my question, physically it means that current goes from right to left and because the off-diagonal is large the current wants to travel along the edges. $\endgroup$
    – Charlie
    Feb 27, 2022 at 2:33

Your Answer

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

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