Help with variational form PDE

Hi, I have a PDE

\begin{align}\frac{\partial P}{\partial L}(L,x,y) &= \frac{A_3^2}{2}\bigg(y^2\frac{\partial^2 P}{\partial x^2}(L,x,y) + x^2\frac{\partial^2 P}{\partial y^2}(L,x,y) + 2xy\frac{\partial^2 P}{\partial x\partial y}(L,x,y)\bigg) \\ &-\bigg(\frac{A_3^2 y^2}{2x} + A_4x\bigg)\frac{\partial P}{\partial x}(L,x,y) - \bigg(\frac{A_3^2 x^2}{2y} - A_4y\bigg)\frac{\partial P}{\partial y}(L,x,y) \\ &+\frac{A_3^2}{2}\bigg(\frac{y^2}{x^2} + \frac{x^2}{y^2} - 2\bigg)P(L,x,y) \end{align}

with initial condition P(L = 0,x,y) = \delta(x-1)\delta(y) and Neumann boundary conditions. The A_3 and A_4 terms are constants.

I rewrite the PDE with \frac{\partial P}{\partial L} \approx \frac{P^{n+1} - P^n}{\Delta L} so

\begin{align} \int_{\Omega}\bigg[v(\mathbf{x})P^{n+1} - \Delta L\bigg[ &\frac{A_3^2}{2}\bigg(y^2\frac{\partial^2 P^{n+1}}{\partial x^2}v(\mathbf{x}) + x^2\frac{\partial^2 P^{n+1}}{\partial y^2}v(\mathbf{x})+ 2xy\frac{\partial^2 P^{n+1}}{\partial x\partial y}v(\mathbf{x})\bigg) \\ &-\bigg(\frac{A_3^2 y^2}{2x} + A_4x\bigg)\frac{\partial P^{n+1}}{\partial x}v(\mathbf{x}) - \bigg(\frac{A_3^2 x^2}{2y} - A_4y\bigg)\frac{\partial P^{n+1}}{\partial y}v(\mathbf{x}) \\ &+\frac{A_3^2}{2}\bigg(\frac{y^2}{x^2} + \frac{x^2}{y^2} - 2\bigg)P^{n+1} v(\mathbf{x})\bigg]\bigg]d\mathbf{x} \\ &= \int_{\Omega}P^nv(\mathbf{x})d\mathbf{x} . \end{align}

where v(\mathbf{x}) is the test function and \mathbf{x} = (x,y). To obtain the weak form I believe I need to remove the second derivatives on P via integration by parts? I do this for the first term by

\int_{\Omega}y^2\frac{\partial^2 P}{\partial x^2}d\mathbf{x} = -\int_{\Omega}y^2\frac{\partial v}{\partial x}\frac{\partial P}{\partial x}d\mathbf{x}

which I think is correct? The second term

\int_{\Omega}x^2\frac{\partial^2 P}{\partial y^2}d\mathbf{x} = -\int_{\Omega}x^2\frac{\partial v}{\partial y}\frac{\partial P}{\partial y}d\mathbf{x}

and I am stuck on the third term

\int_{\Omega}2xyv\frac{\partial^2 P}{\partial x\partial y}d\mathbf{x}

on how to reduce this? I would really appreciate help if what I am doing is correct.

The “correct” integration by parts is partly determined by what is meant by “Neumann boundary conditions”. Dropping boundary terms after integrating by parts implies that certain quantities are zero on the boundary.

I think it would be worth stepping back and looking at whatever physics leads to this PDE. Can it be written in the form of a conservation law,

\frac{\partial P}{\partial L} + \nabla\cdot\mathbf{q} = f\text{ ,}

where \mathbf{q} is some vector-valued flux (depending on P and \nabla P) and f is a source term? (Often the underlying physics gives you this in an integral form, and the PDE comes from divergence theorem and a localization argument.) In conservation laws, a “Neumann boundary condition” usually refers to prescribing the normal component of flux on the boundary. Then the natural way to integrate by parts would be

\int_\Omega v(\nabla\cdot\mathbf{q})\,d\Omega = -\int_\Omega (\nabla v)\cdot\mathbf{q}\,d\Omega + \int_{\partial\Omega}hv\,d\partial\Omega\text{ ,}

where h is the prescribed value of \mathbf{q}\cdot\mathbf{n} from the Neumann boundary condition. (In many cases of practical interest, h is zero, and the boundary integral is dropped, resulting in the “do-nothing” boundary condition, e.g., traction-free boundaries in solid mechanics or insulating boundaries in heat conduction.)

The corresponding conservation law is not immediately clear from the PDE given in your post, but nature doesn’t know about x and y, so there is usually a more compact way to write physically-derived PDEs in terms of \nabla, rather than partials with respect to individual coordinates. (And FEniCS’s UFL then lets you implement solvers directly in that form, using grad, div, etc.)


Hi David, thanks for your reply. I have rewritten the PDE as

\frac{\partial P}{\partial L}(L,x,y) = \nabla\cdot \mathbf{a}^\star\nabla + \mathbf{c}^\star\cdot\nabla P(L,x,y) + QP(L,x,y)


\begin{align} \boldsymbol{a}^\star = \begin{bmatrix} \frac{A_3^2y^2}{2} & \frac{A_3^2xy}{2} \\ \frac{A_3^2xy}{2} & \frac{A_3^2x^2}{2} \end{bmatrix}, \quad \boldsymbol{c}^\star = \begin{bmatrix} -\bigg(\frac{A_3^2y^2}{2x} + A_4x\bigg) - \frac{A_3^2y}{2} \\ -\bigg(\frac{A_3^2x^2}{2y} - A_4y\bigg) - \frac{A_3^2x}{2} \end{bmatrix},\quad Q = \frac{A_3^2}{2}(y^2/x^2 + x^2/y^2 - 2) \end{align}
Approximating the length derivative we have
\begin{align} \frac{P^{n+1} - P^n}{\Delta L} \approx (\nabla\cdot \mathbf{a}^\star\nabla + \mathbf{c}^\star\cdot\nabla)P^{n+1}. \end{align}
Multiply by test function v(\mathbf{x}) and integrate over our domain \Omega to obtain
\begin{align} \int_{\Omega} \bigg(v(\mathbf{x})P^{n+1} - \Delta L\bigg(v(\mathbf{x})\nabla\cdot\mathbf{a}^\star\nabla P^{n+1} + v(\mathbf{x})\mathbf{c}^\star\cdot\nabla P^{n+1} + v(\mathbf{x})QP^{n+1}\bigg)\bigg)d\mathbf{x} = \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}. \end{align}
To formulate the weak form of the problem I integrate by parts the second term to obtain
\begin{align} \int_{\Omega} v(\mathbf{x})\nabla\cdot\mathbf{a}^\star\nabla P d\mathbf{x} = \int_{\Omega} v(\mathbf{x})\mathbf{n}\cdot \mathbf{a}^\star\nabla P ds - \int_{\Omega} \nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P d\mathbf{x}, \end{align}
with Neumann boundary conditions the surface integral vanishes to obtain
\begin{align} \int_{\Omega} v(\mathbf{x})\nabla\cdot\mathbf{a}^\star\nabla P d\mathbf{x} = - \int_{\Omega} \nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P d\mathbf{x}. \end{align}
Hence weak form is written as
\begin{align} \int_{\Omega}\bigg( v(\mathbf{x})P^{n+1} + \Delta L\bigg(\nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P^{n+1} - v(\mathbf{x})\mathbf{c}^\star\cdot\nabla P^{n+1} - v(\mathbf{x})QP^{n+1}\bigg)\bigg)d\mathbf{x} = \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}, \end{align}
where we define
\begin{align} \mathbf{S} = \nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P^{n+1} - v(\mathbf{x})\mathbf{c}^\star\cdot\nabla P^{n+1} - v(\mathbf{x})QP^{n+1}, \end{align}
\begin{align} \int_{\Omega}\bigg( v(\mathbf{x})P^{n+1} + \Delta L\mathbf{S}\bigg)d\mathbf{x} = \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}, \end{align}
\begin{align} a(P,v) = L_{n+1}(v), \end{align}
\begin{align} a(P,v) &= \int_{\Omega}\bigg( v(\mathbf{x})P + \Delta L\mathbf{S}\bigg)d\mathbf{x}, \\ L_{n+1}(v) &= \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}.& \end{align}

My code for solving this is

from __future__ import print_function
from fenics import *
import numpy as np

A_3 = .2
A_4 = .1

#Final L
L = 1
#Number of length steps
num_steps = 10
#Length step size
dl = L / num_steps

#Intial mesh values 
x0 = .1
y0 = 1

# #End mesh values 
x1 = 10
y1 = 10

#Cell numbers in mesh
nx = 100
ny = 100

#Define mesh
mesh = RectangleMesh(Point(x0,y0), Point(x1,y1), nx, ny, "right")

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

#Dirac delta intitial condition
u_D = Expression("(1./(pow(x[0],2) + pow(eps, 2)))*(1./(pow(x[1]-1.,2) + pow(eps, 2)))", eps=1e-6, degree=2)

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

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

#Bilinear form full
a = as_matrix([[.5*A_3**2*x[1]**2, .5*A_3**2*x[0]*x[1]],
[.5*A_3**2*x[0]*x[1], .5*A_3**2*x[0]**2]])

cee = as_vector([-(.5*A_3**2*(x[1]**2/x[0]) + A_4*x[0]) -.5*A_3**2*x[1],
-(.5*A_3**2*(x[0]**2/x[1])) - .5*A_3**2*x[0] ])

#a_grad_u = dot(a,grad(u))
a_grad_u = a*grad(u)
c_dot_grad_u = dot(cee,grad(u))
S = dot(grad(v),a_grad_u) - v*c_dot_grad_u -.5*A_3**2*(x[0]**2/x[1]**2 + x[1]**2/x[0]**2 - 2)*u*v
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 length
    l += dl
    print('L =',l)

    solve(a_form == L_form, u)
    u_file = File('u_sol_.pvd')
    u_file << u, l

    #Probability volume
    print("Area under PDF is:")
    #Compute transmission coefficient at each length step
    print("Transmission coefficient is:")
    print(assemble((u/x[1]**2)*dx ) )

    #Update previous solution
    #Print messages
    print("We are at step",n+1,"out of",num_steps , "step(s)." )

I am interested in computing the quantity in my code assemble((u/x[1]**2)*dx )
but this blows up. Any thoughts? Does my formulation look correct?

So, this derivation assumes that you want a homogeneous Neumann BC on only the “diffusive” flux, q_\text{diff} = -a^\star\nabla P. However, the appropriateness of that depends, again, on the desired physics.

One thing I will point out is that I can’t help but wonder whether some $x$s and $y$s got switched around in the definition of c^\star. The reason I wonder is that

\nabla\cdot\left(c^\star p \right) = c^\star\nabla P + (\nabla\cdot c^\star)P\text{ ,}

and Q is uncannily close to \nabla\cdot c^\star. (If you switch x and y in the A_4 and second A_3 terms from c^\star, it would be equal.) Then you would have a conservation law with

q = -a^\star\nabla P - c^\star P\text{ ,}

where the second term acts like a “advective” flux, in the classic advection–diffusion model. A homogeneous Neumann BC on q\cdot n would preclude all transport of the quantity P through the boundary, and, in the absence of source terms, conserve the integral of P over time.