# Variational form for multi-dimensional diffusion equation

**URL:** <https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606>\
**Category:** variational formulation\
**Created:** [June 22, 2020, 2:43pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606 "2020-06-22T14:43:36Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [June 22, 2020, 2:43pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/1 "2020-06-22T14:43:36Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Yoann](https://avatars.discourse-cdn.com/v4/letter/y/5fc32e/32.png) [@Yoann](https://fenicsproject.discourse.group/u/Yoann)\
**Post date:** [June 22, 2020, 6:57pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/2 "2020-06-22T18:57:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [June 22, 2020, 9:46pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/3 "2020-06-22T21:46:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Yoann](https://avatars.discourse-cdn.com/v4/letter/y/5fc32e/32.png) [@Yoann](https://fenicsproject.discourse.group/u/Yoann)\
**Post date:** [June 23, 2020, 6:31am UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/4 "2020-06-23T06:31:48Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [June 23, 2020, 12:30pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/5 "2020-06-23T12:30:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [June 23, 2020, 2:16pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/6 "2020-06-23T14:16:33Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Yoann](https://avatars.discourse-cdn.com/v4/letter/y/5fc32e/32.png) [@Yoann](https://fenicsproject.discourse.group/u/Yoann)\
**Post date:** [June 23, 2020, 3:08pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/7 "2020-06-23T15:08:25Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [June 23, 2020, 4:21pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/8 "2020-06-23T16:21:02Z")

</div>

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](https://en.wikipedia.org/wiki/Fokker%E2%80%93Planck_equation), 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{ .}

---

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [July 6, 2020, 1:32pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/9 "2020-07-06T13:32:01Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [July 6, 2020, 5:22pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/10 "2020-07-06T17:22:18Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Yoann](https://avatars.discourse-cdn.com/v4/letter/y/5fc32e/32.png) [@Yoann](https://fenicsproject.discourse.group/u/Yoann)\
**Post date:** [July 7, 2020, 9:40am UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/11 "2020-07-07T09:40:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [July 7, 2020, 11:19am UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/12 "2020-07-07T11:19:11Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Yoann](https://avatars.discourse-cdn.com/v4/letter/y/5fc32e/32.png) [@Yoann](https://fenicsproject.discourse.group/u/Yoann)\
**Post date:** [July 7, 2020, 11:30am UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/13 "2020-07-07T11:30:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [July 7, 2020, 12:36pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/14 "2020-07-07T12:36:34Z")

</div>

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}

---

<div class="post-metadata">

**Author:** ![bacsieboy](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bacsieboy/32/1447_2.png) [@bacsieboy](https://fenicsproject.discourse.group/u/bacsieboy)\
**Post date:** [July 7, 2020, 12:51pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/15 "2020-07-07T12:51:57Z")

</div>

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\_3**2_x[2]2,0,A\_32_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])

> ValueError: only 2 non-keyword arguments accepted

Can it be input in matrix form like this?

---

<div class="post-metadata">

**Author:** ![Yoann](https://avatars.discourse-cdn.com/v4/letter/y/5fc32e/32.png) [@Yoann](https://fenicsproject.discourse.group/u/Yoann)\
**Post date:** [July 7, 2020, 1:00pm UTC](https://fenicsproject.discourse.group/t/variational-form-for-multi-dimensional-diffusion-equation/3606/16 "2020-07-07T13:00:22Z")

</div>

> [@bacsieboy](#):
>
> ```auto
> 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]])
> 
> ```

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

```python
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]])

```