How to solve coupled PDEs for electrostatics?


I am looking into solving the 2 following coupled PDEs:
(i) Div[Grad(V)] = -q/cste
(ii) Div[q.Grad(V)] = 0
With q and V being functions: q(x;y) & V(x;y)
With boundary conditions on V all around the 2D domain

  • (i) is a basic Poisson problem and (ii) is a variable coefficient Poisson problem
    Independant resolution of (i) and (ii) works fine with FEniCS, thanks to relevant documentation.
    For both cases: Assuming q is known as a starting point, I can get V as a result on my mesh/domain
  • I tried to couple (i) & (ii) as follows
    STEP1: Assume a first q, set BCs on V and solve (ii) => V1 as an output
    STEP2: Using V1, calculate q with (i) => q2 as an output
    STEP3: Using q2, set (same) BCs on V and solve (ii) => V3 as an output
    Unfortunately, result doesn’t converge
  • I also tried to work with some coupling FEniCS tutorials, but they are more complicated than my need and I couldn’t adapt them

Does anyone know a simple way to couple (i) and (ii) in the same variational statement, and get q and v efficiently ?

In advance Thank you for your help !

I think q = 0 and V given by the poisson equation is a solution.

Guys, do you know how to solve, for example, the potencial of a finite plate? The potencial of the plate is V with regard to the infinity. How can I represent this infinity?

Hi Gabriel,

Your potential decays as V~1/r for a point source. If you look at the largest dimension of your plate, say it is its width W, then basically, it is safe to assume that for any distance r>W, the potential from the plate will decay as V~1/r too (at this scale, your plate is roughly a point source).

What this means is that you do not need an infinite space. You need a domain that is large enough to introduce a minimal mistake on the solution.

Typically you would set up a problem whose analytical solution is known, and you check how large your space needs to be in order to reproduce the solution within some tolerance chosen by you.

1 Like

You are absolutely right.
I just wanted to make sure that there was not a crazy super especial boundary condition.

Thank you for the reply.

Dear Gabriel,

You can find an analytical solution and its numerical implementation in FEniCS for a polarized deformable body, where the far away boundary conditions are set to zero. As a rough estimate, two to three times far from the geometry, it is o.k. to consider “infinity” regarding the boundaries. Have a look at #039 in publications and supply code for the python files under

Best, Emek

Thank you for the trivial solution, but I also need the non trivial one ! :wink:

I made some progresses (at least conceptually) discovering that in order to solve both Poisson Equation and the Transport Equation at the same time, one may use method of characteristics that simplifies the problem. (Charge decays in 1/s along a field line, s being the absicurv). Unfortunately, this is mathematically hard to implement…

If anyone has info about:
how to solve both Poisson equation div(gradV) = q/e0 and Transport equation grad(q.gradV) = 0 in FEniCS, this will be highly appreciated. Thanks !

1 Like

Hello all,
I have been able to progress a bit (at least conceptually) about solving system:
(i) Div[Grad(V)] = -q/cste
(ii) Div[q.Grad(V)] = 0
As we already said, (i) is the basic “hello world” Poisson problem.
(ii) is so called Transport equation, with pre-dominant convection. It may be handled as a variable coefficient Poisson problem, well documented in FEniCS manuals.
Independant solving of (i) and (ii) with FEniCS work fine, but mixing both equations is a convergence nightmare.
1- A conceptual trick is to picture (ii) using the “characteristic method”. When circulating on an E field line, (reminder: E = -Grad(V)), also called a “characteritic curve”, problem simplifies and space charge q decay is given by a simple Ordinary Differential Equation (ODE). Circulating on a E-field line is done using a clever kind of variable change. Has anyone already implemented characteritic method with Python/FEniCS in a mshr grid ?
2- It seems there is also a trick using one of the Galerkin methods. This part is less clear to me… Same question; does anyone have experience on coding this ?
Thank you in advance for sharing experience !

1 Like

“Liked” this post because of its similarity to the problem I am working on currently: setting a far-field boundary condition to be zero in a coupled PDE system returns the trivial result for that variable.

Eagerly watching the space for any help/insight.

Thank You
Warm Regards