Why is the L2 and max error so large when solution is compared to the analytical solution of the heat equation?

variational formulation
1

Hi,
I’ve been trying to validate my solution for an instantaneous point source (IPS) on a finite 2D grid against an analytical solution to an IPS on a 2D infinite domain. I’ve implemented Dirichlet boundary conditions that follows the analytical solution. However, the L2 and maximum error I’m getting is much higher than expected. Any advice would be greatly appreciated.
The analytical solution is shown here: Eq 34.35

#!/usr/bin/env python3
from fenics import *
import numpy as np

#Define sim variables
sim_time = 5.0
time_steps = 50
dt = sim_time/time_steps
mass = 1.0
D_eff = 1.0
x_o = 2.5
y_o = 2.5
tol = 1E-15

#Define mesh and function space
mesh = RectangleMesh(Point(0.0,0.0),Point(5.0,5.0),100,100,'crossed')
V = FunctionSpace(mesh,'P',2)

#Define boundary condition
u_D = Expression('(mass/(4*pi*D_eff*t))*exp((-1/(t*D_eff))*(pow((x[0]-x_o),2)+pow((x[1]-y_o),2)))', \
    degree=5, mass=mass, D_eff=D_eff, x_o=x_o, y_o=y_o, t=1)

def boundary(x,on_boundary):
    return on_boundary

bc = DirichletBC(V,u_D,boundary)

#Initial value of mesh
u_n = interpolate(u_D,V)

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

'''F = u*v*dx + D_eff*dt*dot(grad(u),grad(v))*dx - u_n*v*dx'''
a = (u*v + D_eff*dt*dot(grad(u),grad(v)))*dx
L = u_n*v*dx

#Solving problem
u = Function(V)
t = 1.0
for i in range(time_steps):
    t += dt
    u_D.t = t
    solve(a==L,u,bc)
    error_L2 = errornorm(u_D,u,'L2')
    print("At time: {:.2f} the L2 error is: {:.3f}".format(t,error_L2))
    vertex_values_u_D = u_D.compute_vertex_values(mesh)
    vertex_values_u = u.compute_vertex_values(mesh)
    error_max = np.max(np.abs(vertex_values_u_D-vertex_values_u))
    print("At time: {:.2f} the maximum error is: {:.3f}".format(t,error_max))