Is there an efficient way to use complex Bessel functions in Expression?

Question

I want to use errornorm to compare computed solution with an analytic solution, which uses complex Bessel function of the first kind.

I tried to import them from sympy. If I try to interpolate the Expression over 3D mesh, evaluating takes too much time.

I thought about using Boost, but it doesn't support complex arguments.

Answer 1

Evaluating Python expressions over the entire mesh is costly. Perhaps
precomputations are possible? An example, with Womersley and Bessel
functions can be found in https://bitbucket.org/simula_cbc/cbcflow/
under cbcflow / bcs / Womersley.py.

Answer 2

Hi, as an alternative to Kent's answer consider the following. You might get a meaningful answer from errornorm even if the exact solution is not an Expression. It can be sufficient that it is represented in some higher order space. You can compute expansion coefficients of such a function manually and sympy has tools to speed up the process.

from sympy import symbols, lambdify, besselj
from scipy.special import jv as scipy_besselj
from dolfin import *

mesh = UnitCubeMesh(20, 20, 20)

# Space where the exact solution is represented
Ve = FunctionSpace(mesh, 'DG', 4)
dof_x = Ve.dofmap().tabulate_all_coordinates(mesh).reshape((-1, 3))
X, Y, Z = dof_x[:, 0], dof_x[:, 1], dof_x[:, 2]

# Suppose the solution is besselj(1, x[0]) + besselj(2, x[1]) + besselj(3, x[2])
x, y, z = symbols('x, y, z')
u = besselj(1, x) + besselj(2, y) + besselj(3, z)
u_lambda = lambdify([x, y, z], u, {'besselj': scipy_besselj})
# Interpolate manually
timer = Timer('eval')
timer.start()
u_values = u_lambda(X, Y, Z)
print 'Evaluated %d dofs in % s' % (Ve.dim(), timer.stop())
# Exact solution in Ve
u = Function(Ve)
u.vector()[:] = u_values

# Solution space and the hypot. numerical solution
V = FunctionSpace(mesh, 'CG', 1)
uh = interpolate(Constant(1), V)

print errornorm(u, uh) 

Code above evaluates u_lamba 1.68 million times. On my machine this takes about 3.5 seconds.

Comments

I tried to use this, but I still have a problem with complex numbers:

Traceback (most recent call last):
  File "problem.py", line 9, in <module>
    print u(1.0)
  File "<string>", line 1, in <lambda>
NameError: global name 'I' is not defined

My code (I simplified my formula):

from sympy import I, re, Symbol, lambdify, besselj
from scipy.special import jv 
x = Symbol('x')
u = lambdify(x,(0.00073 - 0.00052*I)*besselj(0,(1.84 - 1.84*I)*x),{"besselj":jv})
# using this line gives the same error:
# u = lambdify(x,(0.00073 - 0.00052j)*besselj(0,(1.84 - 1.84j)*x),{"besselj":jv})
print u(1.0)

---

Try

u = (0.00073 - 0.00052*I)*besselj(0,(1.84 - 1.84*I)*x)
u_lambda = lambdify(x, u, ['numpy', {'besselj': scipy_besselj}]) 
print u_lambda(1.0)

---

Thanks, that does work.