# Using a spline for nonlinear PDE's coefficient

**URL:** <https://fenicsproject.discourse.group/t/using-a-spline-for-nonlinear-pdes-coefficient/6394>\
**Category:** variational formulation\
**Created:** [August 12, 2021, 7:34pm UTC](https://fenicsproject.discourse.group/t/using-a-spline-for-nonlinear-pdes-coefficient/6394 "2021-08-12T19:34:39Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![ankushaggarwal](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ankushaggarwal/32/2846_2.png) [@ankushaggarwal](https://fenicsproject.discourse.group/u/ankushaggarwal)\
**Post date:** [August 12, 2021, 7:34pm UTC](https://fenicsproject.discourse.group/t/using-a-spline-for-nonlinear-pdes-coefficient/6394/1 "2021-08-12T19:34:39Z")

</div>

I want to use a spline function in a PDE. An example is the nonlinear Poisson problem given on FEniCS tutorial [Solving PDEs in Python - \<br\> The FEniCS Tutorial Volume I](https://fenicsproject.org/pub/tutorial/html/._ftut1007.html)  
where instead of the (1+u\*\*2) in q(u), I would like to use a spline (such as a scipy spline interpolation). For example, for a comparison, I was trying

```auto
from scipy.interpolate import UnivariateSpline
import numpy as np
un = np.linspace(0,10,20)
q = UnivariateSpline(un,1+un**2,k=2)

def q_analytical(u):
    "Return nonlinear coefficient"
    return 1 + u**2

#The rest of the code is the same as the tutorial example
from fenics import *
# Use SymPy to compute f from the manufactured solution u
import sympy as sym
x, y = sym.symbols('x[0], x[1]')
u = 1 + x + 2*y
f = - sym.diff(q_analytical(u)*sym.diff(u, x), x) - sym.diff(q_analytical(u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression(u_code, degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = Function(V) # Note: not TrialFunction!
v = TestFunction(V)
f = Expression(f_code, degree=2)
F = q(u)*dot(grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc)

```

However, the above code gives an error because u is a ufl object, while the UnivariateSpline expects a float or numpy array. Any idea on how to solve this?

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [August 13, 2021, 1:16pm UTC](https://fenicsproject.discourse.group/t/using-a-spline-for-nonlinear-pdes-coefficient/6394/2 "2021-08-13T13:16:58Z")

</div>

@kamensky is an expert in this area. But I’ll try not to embarrass myself with the following response:

Splines are awkward to define symbolically since they may be in C^0. Consider the case where you have repeated knots in the knot vector. The `scipy` class you’re using is a numerical function, and obviously shouldn’t be used with symbolic quantities as you’ve pointed out. Just browsing the `sympy` [documentation](https://docs.sympy.org/latest/modules/functions/special.html#b-splines) there’s some support for splines. It might not take too much work to develop your own implementation in `ufl`. Or use `sympy`'s `lamdify` to generate `ufl` expressions.