Variational form for multi-dimensional diffusion equation

variational formulation
1

I’m trying to solve a multi-dimensional diffusion equation (Fokker-Planck equation to be exact) which is multi dimensional (with constants A_{1,2,3})
\frac{\partial P}{\partial L}(L,x,y,z) = \frac{A_3^2z^2}{2}\frac{\partial^2 P}{\partial x^2}(L,x,y,z) \\ + A_3^2xz\frac{\partial^2 P}{\partial x\partial z}(L,x,y,z) + \bigg(A_1^2+A_2^2+A_3^2\frac{z^2}{x^2}\bigg)\frac{\partial^2 P}{\partial y^2}(L,x,y,z) \\ + A_3^2x^2\frac{\partial^2 P}{\partial z^2}(L,x,y,z) +\bigg(\frac{A_3^2}{2}-A_4\bigg)x\frac{\partial P}{\partial x}(L,x,y,z)\\ + \bigg(\frac{A_3^2}{2} +A_4\bigg)z\frac{\partial P}{\partial z}(L,x,y,z) - A_5\frac{\partial P}{\partial y}(L,x,y,z) \\ = F(L,x,y,z).

PDE has initial condition P(L=0,x,y,z) = \delta(x)\delta(y)\delta(z-1) with x,\,y\in[1,\infty) z \in [1,\infty) and Neumann as standard. The L variable can be thought of as time (it corresponds to length). First I discretise the left hand side via
\bigg(\frac{\partial P}{\partial L}\bigg)^{n+1} \approx \frac{P^{n+1}-P^n}{\Delta L}

to obtain

P^{n+1} - P^n = \Delta L \,F(L,x,y,z).

My question is now how to obtain the correct variational form. What should my test function depend on? Should it be multidimensional too? I am thinking I should have something like

\int_{\Omega}\bigg(P(L,x,y,z)v(\mathbf{x}) - \Delta L\, F(L,x,y,z)v(\mathbf{x})\bigg)\,d\mathbf{x} = \int_{\Omega}p^{n}(L,x,y,z)v(\mathbf{x})\,d\mathbf{x}

where \mathbf{x} = [x,y,z]^{T}, but I am really confused on how to reduce the second derivatives in this case. Thank you for help in advance.

2

First do you have any boundary conditions ? It will help to reduce the bilinear form. If you have Dirichlet boundary conditions then you can proceed to the reduction of the second derivatives . For example:
\int \frac{A_3^2z^2}{2}\frac{\partial^2 P}{\partial x^2}(L,x,y,z) v(\mathbf{x})\mathrm{d} \mathbf{x}= \int \frac{A_3^2z^2}{2}\frac{\partial P}{\partial x}(L,x,y,z) v(\mathbf{x})n_x\mathrm{d} S -\int \frac{\partial P}{\partial x}(L,x,y,z) \frac{\partial}{\partial x}\left(\frac{A_3^2z^2}{2}v(\mathbf{x})\right)\mathrm{d} \mathbf{x}
Please take care that \mathbf{n} is oriented towards the outside of the surface.

Your space should be \mathbb{R}^3 for the variables (x,y,z) and the discretisations that you proceeded for the dimension of L is similar to time discretisation in the tutorials. Thus to solve your PDE you have to do a \texttt{while} loop over L.
Your test function and trial function should be declare on \mathbb{R}^3

3

Hmm. No Dirichlet BC’s, only Neumann at the boundary going to zero. I’m not sure the form for the test function is correct. I’m still unsure if this will work.

4

When you say

I’m not sure the form for the test function is correct.

Are you talking about the right hand side ? If it’s the case, yes your form is correct

L(v)=\int_\Omega P^n(L,x,y,z)v(\mathbf{x})\,\mathrm{d}\mathbf{x}

If you apply Neuman boundary conditions then you have :

\int_\Omega \frac{A_3^2z^2}{2}\frac{\partial^2 P}{\partial x^2}(L,x,y,z) v(\mathbf{x})\mathrm{d} \mathbf{x}= -\int_\Omega \frac{\partial P}{\partial x}(L,x,y,z) \frac{\partial}{\partial x}\left(\frac{A_3^2z^2}{2}v(\mathbf{x})\right)\mathrm{d} \mathbf{x}

And you can do the same for each term since \nabla P =0 on \partial \Omega.
I think it can work because you’re applying initial conditions. Otherwise, if it was a simple PDE then it would be ill posed because P and P+C are solution… Don’t you have a conservation equation or something like this?

5

Hi Yoann, merci beacoup. I think this may be correct actually. I will try this.

6

Hi Yoann, small question. Does the \int_{\Omega} A_3^2xz\frac{\partial^2 P}{\partial x\partial z}v(\mathbf x)\,d\mathbf x term need to be reduced? I’m not sure how to reduce this term since it contains derivative in x and z?

7

De rien, (Are you french either by the way? )

Yes you have to reduce this term. If you suppose that P is \mathscr{C}^2 then you can choose whether x or z thanks to Cauchy-Schwarz, otherwise you should do the integration by part over x.

8

Can the strong form of the PDE be re-written in the form of the partial time derivative equal to the divergence of some flux? In that case, the obvious “correct” choice for integration by parts would be to integrate by parts to make the divergence into a gradient on the test function. I know nothing about the Fokker–Planck equation, but, in the Wikipedia article, under “Higher dimensions”, both forms listed have a spatial divergence operator as the “outer-most” operation on the right-hand side, so you could do

\cdots = \int_\Omega \sum_i\left\lbrack\frac{\partial}{\partial x_i}(\cdots)_i\right\rbrack v\,d\Omega = -\int_\Omega\sum_i\left\lbrack(\cdots)_i\frac{\partial v}{\partial x_i}\right\rbrack\,d\Omega + (\text{bdry. terms}) \text{ .}
9

I can write the PDE as

\frac{\partial P}{\partial t}(L,x_1,x_2,x_3) = \bigg(\sum_{i,j=1}^{3}a_{i,j}\frac{\partial^2}{\partial x_i\partial x_j} + \sum_{i=1}^{3}b_{i}\frac{\partial}{\partial x_i}\bigg)P(t,x_1,x_2,x_3)

where

a_{ij} = \begin{bmatrix} A_3^2x_3^2 & 0 & A_3^2x_1x_3 \\ 0 & (A_1^2+A_2^2+A_3^2\frac{x_3^2}{x_1^2}) & 0 \\ A_3^2x_1x_3 & 0 & A_3^2x_1^2 \end{bmatrix}
b = \begin{bmatrix}A_4x_1 \\ A_5 \\ \frac{A_3^2}{2}x_3 - A_4x_3 \end{bmatrix}

so I don’t think it’s possible?

10

I don’t immediately see a way to reformulate the right-hand side in terms of a flux divergence, but it may be worth digging further back into how this equation was derived. In particular, without a conservation law structure to the governing PDE, do you have any other mathematical argument (aside from P nominally being a probability distribution) for why \int P\,d\mathbf{x} should be constant?

1 Like
11

You can write the PDE as:
\partial_t P=(A\nabla+b)\cdot \nabla P
But to simplify A\nabla+b you need relations between the coefficients A_i.
Still you can write (A\nabla+b)\cdot \nabla =\nabla\cdot R\nabla + c\cdot \nabla
R=\left[\begin{array}[ccc] eA_3^2x_3^2 & 0&A_3^2x_3x_1 \\ 0&A_1^2+A_2^2+A_3^2\frac{x_3^2}{x_1^2}&0\\ A_3^2x_3x_1 & 0&A_3^2x_1^2 \end{array}\right]
c=\left[\begin{array}[c] eA_4x_1-A_3^2x_3\\ A_5\\ -\frac{A_3^2}{2}x_3-A_4x_3 \end{array}\right]
You must check but I think it’s correct.
Then what you have is:
\int v \, \nabla\cdot R\nabla P\, \mathrm{ d}x=\int v \, \mathbf{n}\cdot R\nabla P \,\mathrm{ d}s-\int \nabla v\cdot R\nabla P \,\mathrm{ d}x
Then with BCs
\int v \, \nabla\cdot R\nabla P\, \mathrm{ d}x=-\int \nabla v\cdot R\nabla P \,\mathrm{ d}x

1 Like
12

Hi Yoann, thanks for your reply. What is A? Is it my a_{i,j} matrix? Also I’m not sure where the -A_3^2z term comes from in the first entry of c and also where the negative sign change from the last term in c comes from?

13

yes A is your “a_{i,j}” matrix. The terms are just a way to rewrite your PDE what I’ve done is simply:
\sum_{i,j}a_{i,j}\partial_i \partial_j +\sum_{i}b_{i}\partial_i =\sum_{i,j}\partial_i r_{i,j} \partial_j +\sum_{i}c_{i}\partial_i

c=b+\left[\begin{array}[c] 0-A_3^2x_3\\ 0\\ -A_3^2x_3 \end{array}\right]
You can find c by studying each term of A.

1 Like
14

Hi. I tried this. I obtained

\int_{\Omega}vp^{n+1} + \Delta L\bigg(\nabla v\cdot \underbrace{R\nabla p^{n+1}}_{R\,dot\,grad\,u} - v\underbrace{c\cdot\nabla p^{n+1}}_{c\,dot\,grad\,u}\bigg)d\mathbf{x} = \int_{\Omega}p^nvd\mathbf{x}

15

I tried implementing

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt


L = 2   #Final L
num_steps = 5  #Number of length steps
dl = L / num_steps #Length step size

#Define positive constants
A_1 = 0.0001
A_2 = 0.0001
A_3 = 0.0001
A_4 = 0.0001
A_5 = 0.0001

#Intial mesh values
x0 = 0.01
y0 = 0.01
z0 = 1

#End mesh values
x1 = 2
y1 = 2
z1 = 2

#Number of cells in each direction of the mesh
nxyz = 10


#Define mesh
mesh = BoxMesh(Point(x0, y0, z0),Point( x1, y1, z1), nxyz, nxyz, nxyz)

V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
mpx=np.array(mesh.coordinates())
dx = Measure('dx',mesh)

#Dirac Delta (mass) initial condition with epsilon definition
u_D = Expression("(1./(pow(x[0],2) + pow(eps, 2)))*(1./(pow(x[1],2) + pow(eps, 2)))*(1./(pow(x[2]-1.,2) + pow(eps, 2)))", eps=1e-6, degree=2)


#Define initial value
u_n = interpolate(u_D, V)
C = assemble(u_n*dx)
u_n.assign(u_n/C)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

#Bilinear form
R = np.array([A_3**2*x[2]**2,0,A_3**2*x[0]*x[2]],[0,(A_1**2+A_2**2+A_3**2*(x[2]/x[0])**2),0],[A_3**2*x[0]*x[2],0,A_3**2*x[0]**2])
c = np.array([-A_3**2*x[2]+A_4*x[0]], [A_5], [0.5*A_3**2*x[2]-A_4*x[2]])

R_dot_grad_u = R.dot(grad(u))
c_dot_grad_u = c.dot(grad(u))
S = dot(grad(v),R_dot_grad_u) - dot(v,c_dot_grad_u)
a_form = (u*v + dl*(S))*dx
L_form = u_n*v*dx

#Starting to build length steps
u = Function(V)
l = 0


for n in range(num_steps):

    #Update current lengthj
    l += dl

    #Solve
    solve(a_form == L_form, u)

but there is an error > R = np.array([A_32x[2]2,0,A_32x[0]*x[2]],[0,(A_12+A_22+A_32*(x[2]/x[0])2),0],[A_32*x[0]x[2],0,A_3**2x[0]**2])

ValueError: only 2 non-keyword arguments accepted

Can it be input in matrix form like this?

16

Replace this by the following to correctly declare your numpy array:

R = np.array([[A_3**2*x[2]**2,0,A_3**2*x[0]*x[2]],[0,(A_1**2+A_2**2+A_3**2*(x[2]/x[0])**2),0],[A_3**2*x[0]*x[2],0,A_3**2*x[0]**2]])
c = np.array([-A_3**2*x[2]+A_4*x[0], A_5, 0.5*A_3**2*x[2]-A_4*x[2]])
1 Like