Solving a pde
\frac{\partial p}{\partial L}(L, \eta)= \frac{\partial}{\partial \eta}\left[\left(\eta^{2}-1\right) \frac{\partial p}{\partial \eta}(L, \eta)\right], \quad \eta>1,
with Dirac delta initial condition
p(L=0, \eta)=\delta(\eta-1),
and Dirichlet bc P(L,\eta=1) = 0.
Variational form
a(P, v)=\int_{\Omega}\left(P v+\frac{\Delta L}{\theta}\left(\eta^{2}-1\right) \frac{\partial P}{\partial \eta} \frac{\partial v}{\partial \eta}\right) d \eta,
L(P,v) = \int_{\Omega}Pv\,d\eta.
Code for this
"""
1d Fokker-Planck equation, spatially dependant diffusion process with
Dirac delta initial condition, Dirichlet and Neumann.
"""
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
L = 100 # final L
num_steps = 5000 # number of length steps
dl = L / num_steps # length step size
theta = 1 #Lloc constant term
# Create mesh and define function space
nx = 100
mesh = IntervalMesh(nx, 1., 1000)
#mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
# Define boundary condition
u_D = Expression("1./(pow(x[0]-y,2) + pow(eps, 2))", eps=1e-6, pi=pi, y=1., degree=2)
#BC Using tanh approximation
#u_D = Expression("7./pow(cosh(7.*(x[0]-1.)),2)", eps=1e-6, pi=pi, y=1., degree=2)
# plot(u_D)
def right(x, on_boundary):
return near(x[0], 1.)
def left(x, on_boundary):
return near(x[0], 1.)
bc = DirichletBC(V, Constant(0), right)
# 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)
a_form = (u*v + (dl/theta)*(x[0]*x[0]-1)*u.dx(0)*v.dx(0))*dx
L_form = u_n*v*dx
# Time-stepping
u = Function(V)
l = 0
#List for values of u(eta)
u_eta = []
L_vals = []
mean_t = []
for n in range(num_steps):
# Update current time
l += dl
# Compute solution
solve(a_form == L_form, u, bc)
vtkfile = File('foquesol.pvd')
vtkfile << u
#Compute moments
mean_t.append(assemble((2/(1+x[0]))*u*dx))
print(assemble((2/(1+x[0]))*u*dx))
#Plot solution
#plot(u)
# Update previous solution
u_n.assign(u)
#print(assemble(u*dx))
#Store u(eta) values into list and also corresponding lengths
#u_eta.append(u(l))
L_vals.append(l)
#print(l)
# #plot PDF solution
plt.scatter(L_vals,mean_t)
plt.xlabel('$L$')
plt.ylabel('$\\mathbb{E}(\\tau)$')
plt.show()
# Hold plot
#print(len(mean_t))
Once I solve for P(L,\eta) I compute moments of this function, since P is actually a probability density function. The moment integral is I = \int_{1}^{\eta_{max}}\frac{2}{1+\eta}P(L,\eta)\,d\eta which is calculated at each l inside the loop with mean_t. I have the analytical result for this pde which is in the form of an integral solution. Plotting the analytical solution gives
My problem is that the plot I get of I(L) in FeniCS looks like this
I should get I(L=0) = 1 like in the analytical result. I think there is a problem implementing the initial Dirac mass. I have tried using epsilon delta method provided by Jeremy B and also via
\delta(\eta-1) \approx \frac{n}{\cosh{(n(\eta-1))}^2}, \quad n >,1
with n=7 in my code above. Both give the same problem. If anyone could provide me with some advice or things to try to get these solutions to match up that would be great.