I’m solving a PDE with Dirac delta initial condition
\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{L_{loc}}\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). With the help of Jeremy B I have code which solves 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
L = 2 # final L
num_steps = 10 # 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., 3)
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)
# 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
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
# Plot solution
plot(u)
#plot(mesh)
# Update previous solution
u_n.assign(u)
print(assemble(u*dx))
# Hold plot
plt.xlabel('$\\eta$')
plt.ylabel('$P(\\eta,L)$')
plt.legend(['$P(\\eta,L_0)$','$P(\\eta,L_1)$'], loc='upper right')
plt.show()
I want to use the solution from this code to calculate statistical moments, namely the integral
\int_{1}^{\infty}\left(\frac{2}{1+\eta}\right)^{n} P(L, \eta) d \eta
I know how to do this in python using scipy where I know what the function P is, but not in this case where it is a vector. The resulting plot I want is something like this, where I plot the above integral with varying L on the x axis
Sorry if this is a really dumb question but I’m stuck on the best way to try and achieve this.