 # Curve_fit of exponential function using FEM solution as input

Hi everyone,

I have fit an exponential function to some data, and now I want to use this exponential function with T, the solution to my finite element simulation, as an input. However, I am encountering the following error:

``````File "target-spotsize-fixedP.py", line 47, in exponenial_func
return a*exp(-b/x)
ufl.log.UFLValueError: Invalid type conversion: [...] can not be converted to any UFL type

``````

However I cannot see where my problem is, is there something wrong with using Project() in this case to output the value of this exponential function across my whole mesh?

Below is my simplified code:

``````from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from ufl import as_tensor
from ufl import Index
import math
from scipy.optimize import curve_fit
from mpi4py import MPI

mesh = UnitSquareMesh(16,16)

x = [0.01,10,20,30,40,50,60,80,100,1000,10000]  #T range
y = [0,0.1,0.36,0.61,0.72,0.78,0.82,0.84,0.87,0.99,0.99999] # y to fit

def exponenial_func(x, a, b):
return a*np.exp(-b/x)

popt1, pcov1 = curve_fit(exponenial_func, x, y ,p0 = [2, 10])

Space = FunctionSpace(mesh,'P', 1)
T = project(Expression("x+0.1*x", degree=2), Space)

vp_space = FunctionSpace(mesh, 'CG',1)
vp_out = Function(vp_space, name = '\$T\$')
vp_ = exponenial_func(T,*popt1)
vp_out = project(vp_, vp_space,solver_type='cg')

``````

`exp` is overloaded (where possible, in general never use `from some_module import *`).
Since `curve_fit` belongs to scipy, you want to use `scipy.exp` (or related, such as `numpy.exp`). At the moment your code is using `ufl.exp`. scipy doesn’t know about ufl.

Unfortunately using `scipy.exp` or `np.exp` gives me another error, namely:

`````` File "target-spotsize-fixedP.py", line 243, in <module>
vp_ = exponenial_func(T,*popt1)
File "target-spotsize-fixedP.py", line 48, in exponenial_func
return a*scipy.exp(-b/x)
TypeError: loop of ufunc does not support argument 0 of type Division which has no callable exp method
``````

You’ll need to provide a minimal working example of the problem. Your code does not work in any case,

``````NameError: name 'temp_vp' is not defined
``````

Ah, thanks for pointing that out. I have altered the MWE accordingly.

It’s not really clear what you’re trying to do. You’re setting `vp_ = exponenial_func(T,*popt1)`, but `exponenial_func` is a normal python function (taking “x” as first argument) and `T` is declared as an undefined dolfin Function, `T = Function(Space)`. What could `exponenial_func(T)` possibly mean?

What I am trying to do here is take the solution T, from across my mesh and then calculate the value of `exponenial_func` according to the fit calculated earlier with the solution, T at a point on the mesh, taken as the x argument in the function. My hope would then this would calculate the value of `exponenial_func` across the whole mesh.

I think I see what you mean. But your code is not doing that. You are not providing any earlier solution to your `exponenial_func`. You’ve declared `T` as an empty (undefined) function, `T = Function(Space)` . It has no value.

Also your sample code does not run anyway (never gets to the line you’re trying to fix). Please make sure you’ve run and tested your code.

Ok, so I have updated the code slightly, now it is at the line that I wish to fix. I have also altered `exp()` to `np.exp()` but there is still this issue that I quoted in the original question.

Right, so now your dolfin function `T` is defined, and you want to pass it into your `exponenial_func`.

If I understand correctly, your situation is that you have some parametric function f(t) (your exponential function), and some other function T(x) (your dolfin function). What you’re trying to do is generate a third function in x,

g(x) = f( T(x) )

The title in your initial question was confusing, it sounded like you want to apply curve fitting to your calculated solution T(x).

Probably what you want to do is set up your transform function as an Expression, not a python function. Then you can project it like you did with your `T` function. But the Expression needs to take T as a parameter. That’s no longer simple with a simple in-line `Expression`, instead try a `UserExpression`

``````class exponenial_expr(UserExpression):
def __init__(self, t, a, b, **kwargs):
super().__init__(**kwargs)
self.t = t
self.a = a
self.b = b
def eval(self, values, x):
values = exponenial_func( self.t(x), self.a, self.b )
def value_shape(self):
return ()
``````

Then you can project,

``````vp = project( exponenial_expr( T,*popt1 , degree=2), Space)
``````

or you might prefer to interpolate

``````vp_interp = interpolate( exponenial_expr( T,*popt1 , degree=2), Space)
``````
1 Like

Wonderful! That works excellently, thank you very much!