Formulation weak form of PDE for solving in FeniCS I'm trying to derive the weak form for the following Fokker-Planck equation to use in a finite element package. I have
$$\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{\theta}\bigg(2\eta\frac{\partial P}{\partial \eta} +(\eta^2+1)\frac{\partial^2 P}{\partial \eta^2}\bigg), \quad \eta >1,$$
With initial condition $P(L=0) = \delta(\eta-1)$. The Dirichlet boundary conditions can be set later. I require the problem to be in the weak form (page 11) for it to be solved by the FeniCS finite element package. I choose to discretize the time derivative via
$$\bigg(\frac{\partial P}{\partial L}\bigg)^{n+1} \approx \frac{P^{n+1}-p^n}{\Delta L}$$
so
$$P^{n+1}-\frac{\Delta L}{\theta}\bigg(2\eta\frac{\partial P^{n+1}}{\partial\eta}+(\eta^2+1)\frac{\partial^2 P^{n+1}}{\partial \eta^2}\bigg) - p^n = 0.$$
The weak form of the problem I want to write like
$$a(P,v) = L_{n+1}(v)$$
Where $v$ is a test function. I think then I can write
$$\begin{align}
a(P,v) &= \int_{\Omega}P^{n+1}v-\frac{\Delta L}{\theta}\bigg(2\eta\frac{\partial P^{n+1}}{\partial\eta}+(\eta^2+1)\frac{\partial^2 P^{n+1}}{\partial \eta^2}\bigg)\,dx, \\
L_{n+1}(v) &= \int_{\Omega}p^{n}v\,dx.
\end{align}$$
I'm not sure whether this is correct, and if I should do something with the second derivative in $a(P,v).$
 A: Actually, the boundary conditions may change the expression of the weak form. Here I am going to assume that you have homogeneous Dirichlet boundary conditions at $\eta=1$ as well as at $\eta_\max>1$ which is where your are going to truncate your (at the moment unbounded, so untractable numerically) domain in $\eta$.
Observe that your equation rewrites
$$
\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{\theta}\frac{\partial}{\partial \eta}\bigg( (\eta^2+1)\frac{\partial P}{\partial \eta}\bigg)
$$
So after your backward Euler discretization in time, it yields
$$
P^{n+1} - \frac{\Delta L}{\theta}\frac{\partial}{\partial \eta}\bigg( (\eta^2+1)\frac{\partial P^{n+1}}{\partial \eta}\bigg) = P^n
$$
in strong form.
To obtain the weak form you can formally integrate the equation against a test function $v$ that vanishes at the endpoints of your interval. It yields after one integration by parts in the second term of the left-hand side
$$
\int_1^{\eta_\max}P^{n+1}(\eta)v(\eta)d\eta + \int_1^{\eta_\max}\frac{\Delta L}{\theta}(\eta^2+1)\frac{\partial P^{n+1}}{\partial \eta}(\eta)\frac{\partial v}{\partial  \eta}(\eta)d\eta = \int_1^{\eta_\max}P^n(\eta)v(\eta)d\eta
$$
From there you have your bilinear form
$$
a(P,v)=\int_1^{\eta_\max}P(\eta)v(\eta)d\eta + \int_1^{\eta_\max}\frac{\Delta L}{\theta}(\eta^2+1)\frac{\partial P}{\partial \eta}(\eta)\frac{\partial v}{\partial  \eta}(\eta)d\eta
$$
and your linear form
$$
L_n(v)=\int_1^{\eta_\max}P^n(\eta)v(\eta)d\eta
$$
