# Software tutorial/Functions as objects

# Background

From this point onward in the course we will often have to deal with an arbitrary function, such as finding the function's zeros, or integrating the function between two points. It is helpful if we can write MATLAB or Python code that can operate on **any** function, not just a specific function, \(f(x)\).

For example, if we write a Python function to find the zero of \(f(x) = 3x - \frac{2}{x}\), then we would like to send that function, let's call it `my_function`, into another Python function that will find the zero of any function, \(f(x)\).

```
def my_function(x):
"""
Returns the value of f(x) = 3x - 2/x, at the given x
"""
return 3*x - 2/x
lower = 0.1
upper = 3.0
root = bisection(my_function, lower, upper)
```

Python (or MATLAB) will see that `my_function` isn't an ordinary variable, it is another function that will return the value of \(f(x)\) when given a value \(x\). This means that when the code inside the `bisection` routine wants to evaluate \(f(x)\), it can do so. If you want to change which function is being operated on, then you just call `bisection`, but change the first input to point to a different function.

## MATLAB details: inline functions, and function handles

MATLAB has 3 ways to pass a function into another function. For simple, one-line functions, you can use inline or anonymous functions. For more complex functions you will use function handles.

### Inline functions

For a simple, one-line function, it is easy to create and use an *inline* function:

```
>> f = inline('3*x - 2./x')
f =
Inline function:
f(x) = 3*x - 2./x
>> f(1.5)
ans =
3.1667
>> f([0, 1.5, 3.0])
ans =
-Inf 3.1667 8.3333
```

The variable `f` is an object, of the type `inline`:

```
>> whos f
Name Size Bytes Class Attributes
f 1x1 836 inline
```

You can construct more complex inline functions. For example, here we construct a function of two variables, \(f(x,y) = 3x^2 - 4x^3 - 2xy\)

```
>> f = inline('3.*x.^2 - 4.*x.^3 - 2.*x.*y', 'x', 'y') % tell MATLAB what the variable names are
f =
Inline function:
f(x,y) =
>> f(3, 4)
ans =
-105
```

### Anonymous functions

These functions are very similar to inline functions. Continuing the example above:

```
>> f = @(x, y) 3.*x.^2 - 4.*x.^3 - 2.*x.*y
f =
@(x,y)3.*x.^2-4.*x.^3-2.*x.*y
>> f(3,4)
ans =
-105
>> f([3, 2, 3, 2], [4, 4, 5, 5])
ans =
-105 -36 -111 -40
```

The first instruction creates an anonymous function, called `f`, which accepts two inputs, `x` and `y`. You can now call this function, `f`, and give it either scalar inputs, or even vector inputs. In the last case you will get vector outputs.

The syntax for anonymous functions is always: `@(`*one or more comma-separated variable names*`)` *a valid MATLAB expression*.

### Function handles

Later on in the course you will need to construct more complex functions that cannot be written on a single line. For example, we derived the heat balance for a jacketed CSTR in which a reaction takes place.

```
% This code appears in a file called "jacketed_CSTR.m"
function dTdt = jacketed_CSTR(T, CA, Tin, Tj)
% Calculates the time-derivative of temperature in a jacketed CSTR.
%
% INPUTS
% T : tank temperature [K]
% CA : concentration of A in tank [mol/m^3]
% Tin : inlet temperature [K]
% Tj : jacket temperature [K]
% OUTPUT
% dTdt: time-derivative of tank temperature [K/s]
vol = 5.5; % m^3 Tank volume
U = 6000; % W/(m^2.K) Heat transfer coefficient
A = 15; % m^2 Jacket surface area
Fin = 0.25; % m^3/s Inlet flowrate
rho = 1140; % kg/m^3 Liquid phase density
Cp = 4300; % J/(kg.k) Liquid heat capacity (assumed constant)
k0 = 0.004; % mol/(m^3.s) Reaction rate
Hr = 65800; % J/mol Heat of reaction
E = 20000; % J/mol Activation energy
R = 8.314; % J/(mol.K) Gas constant
enthalpy = Fin/vol .* (Tin - T);
jacket = -U * A .* (T - Tj) ./ (rho * Cp * vol);
reaction = k0*(-Hr)/(rho*Cp) .* CA .* exp(-E./(R.*T));
dTdt = enthalpy + jacket + reaction;
```

Now we would like to use the above function, which returns the time-based derivative, \( \displaystyle \frac{dT(t)}{dt} \), and find the tank temperature, \(T\) that sets this derivative is zero.

But, we can see the function above has 4 inputs, but we are only going to vary the first input, \(T\). How do we specify the other 3 inputs? We will create an anonymous function, using a function handle:

```
>> concA = 0.4; % mol/m^3
>> Tin = 300; % K
>> Tj = 350; % K
>> heat_balance = @(T)jacketed_CSTR(T, concA, Tin, Tj);
>> whos
Name Size Bytes Class Attributes
Tin 1x1 8 double
Tj 1x1 8 double
concA 1x1 8 double
heat_balance 1x1 16 function_handle
```

So we have fixed the other 3 constants to the given values, and created an anonymous function using the `@(...)` syntax. However, in this case, the part after the `@(...)` is not a MATLAB expression, it is the name of the function file we wish to call. In this case, we wish to call `jacketed_CSTR.m`, with the four inputs.

Let's use this function handle now:

```
>> T = [290, 300, 310, 320, 330];
>> heat_balance(T)
ans =
0.6548 0.1669 -0.3210 -0.8089 -1.2969
```

This shows that the function's output changes sign somewhere between \(T=300\)K and \(T=310\)K. We can even use function handles in plot commands. Try this:

```
>> T = 300 : 0.1 : 310;
>> plot(T, heat_balance(T))
>> grid on
```

and use the graph to locate the function's zero.

## Python details: functions are objects

Things are a little simpler in Python: there is only one type of function. But to understand that you must also understand the concept of an `object`. When you write `a = 42` in Python, it creates an object, of type integer, and gives it the name, called `a`.

**Simple rule**: *Everything is an object in Python*.

All four of these variables, `my_string, my_integer, my_float` and `my_function` are objects.

```
my_string = 'abc'
my_integer = 123
my_float = 45.6789
def my_function(x):
return 3*x - 2/x
```

You can use the `type(...)` function to determine an object's type:

```
>>> type(my_string) # string object
<type 'str'>
>>> type(my_integer) # integer object
<type 'int'>
>>> type(my_float) # floating point object
<type 'float'>
>>> type(my_function) # function object
<type 'function'>
```

Python has many built-in types, and you can even create your own. Let's repeat the MATLAB example above for the jacketed CSTR:

```
# Copy and paste this code into a .py file and run it.
import numpy as np
def jacketed_CSTR(T, CA, Tin, Tj):
"""
Calculates the time-derivative of temperature in a jacketed CSTR.
INPUTS
T : tank temperature [K]
CA : concentration of A in tank [mol/m^3]
Tin : inlet temperature [K]
Tj : jacket temperature [K]
OUTPUT
dTdt: time-derivative of tank temperature [K/s]
"""
vol = 5.5 # m^3 Tank volume
U = 6000 # W/(m^2.K) Heat transfer coefficient
A = 15 # m^2 Jacket surface area
Fin = 0.25 # m^3/s Inlet flowrate
rho = 1140 # kg/m^3 Liquid phase density
Cp = 4300 # J/(kg.k) Liquid heat capacity (assumed constant)
k0 = 0.004 # mol/(m^3.s) Reaction rate
Hr = 65800 # J/mol Heat of reaction
E = 20000 # J/mol Activation energy
R = 8.314 # J/(mol.K) Gas constant
enthalpy = Fin/vol * (Tin - T)
jacket = -U * A * (T - Tj) / (rho * Cp * vol)
reaction = k0*(-Hr)/(rho*Cp) * CA * np.exp(-E/(R*T))
return enthalpy + jacket + reaction
if __name__ == '__main__':
print(jacketed_CSTR) # <function jacketed_CSTR at 0x42d1f0>
print(type(jacketed_CSTR)) # <type 'function'>
concA = 0.4 # mol/m^3
Tin = 300 # K
Tj = 350 # K
T = np.array([290, 300, 310, 320, 330])
output = jacketed_CSTR(T, CA=concA, Tin=Tin, Tj=Tj)
print(output) # [ 0.6548 0.1669 -0.3210 -0.8089 -1.2969]
```

# Find zeros of a nonlinear function

So at this stage we assume you have a function file that you want to find the zero for. For MATLAB users you must have a function handle, or for Python users, you must have just a regular function.

MATLAB users only have one built-in routine to find the zero of a function: `fzero(...)`. This MATLAB function starts with an open method (it only requires a single starting point), but will automatically switch to the bisection method if the open methods fail. Newton's method, the secant method, and the bisection method are **not available** as built-in MATLAB functions: you must write your own code if you need this.

Python users have several functions available, depending on the problem. The methods we have learned about in this course are all implemented already:

`bisect`: implements the bisection method.`fixed_point`: implements the fixed-point algorithm`brentq`: implements Brent's method. This is probably the closest to MATLAB's`fzero`function, but without the automatic bisection method.`newton`: implements both Newton's method and the Secant method (see below)

We will first present the MATLAB code, then the Python code.

MATLAB |
---|

```
% Continuing our example from above:
>> concA = 0.4; % mol/m^3
>> Tin = 300; % K
>> Tj = 350; % K
>> heat_balance = @(T)jacketed_CSTR(T, concA, Tin, Tj);
% Find the zero, starting from an initial guess of 300K
>> fzero(heat_balance, 300)
ans =
303.4208
% If you like more information about how it found the zero
>> fzero(heat_balance, 300, optimset('Display','iter'))
Search for an interval around 300 containing a sign change:
Func-count a f(a) b f(b) Procedure
1 300 0.166908 300 0.166908 initial interval
3 291.515 0.580928 308.485 -0.247112 search
Search for a zero in the interval [291.515, 308.485]:
Func-count x f(x) Procedure
3 308.485 -0.247112 initial
4 303.421 1.15945e-10 interpolation
5 303.421 -5.25757e-16 interpolation
6 303.421 -5.25757e-16 interpolation
Zero found in the interval [291.515, 308.485]
ans =
303.4208
``` Shows that it constructed an interval of \([291.515, 308.485]\), and then used an interpolation method to find the root, \(T=303.4\)K. |

Python |
---|

```
>>> import numpy as np
>>> from scipy.optimize import *
>>> concA = 0.4 # mol/m^3
>>> Tin = 300 # K
>>> Tj = 350 # K
# Apply the bisection method:
>>> bisect(jacketed_CSTR, 300, 310, args=(concA, Tin, Tj))
303.42075240691372
# Apply the Secant method: use the ``newton`` function, but don't give it f'(x)
>>> newton(jacketed_CSTR, 300, args=(concA, Tin, Tj))
303.42075240691321
# Apply Brent's method (similar to fzero in MATLAB):
>>> brentq(jacketed_CSTR, 300, 310, args=(concA, Tin, Tj))
303.42075240691321
``` In all the cases above you can get a bit more information about how the zero was found, just add an extra input to each function: ```
>>> root, solution = bisect(jacketed_CSTR, 300, 310, args=(concA, Tin, Tj), full_output=True)
>>> solution
<scipy.optimize.zeros.RootResults object at 0x45d3bb0>
# So we get a new object as the output, a ``RootResults`` object. You can inspect this object:
>>> solution.converged
True
>>> solution.root
303.42075240691372
>>> solution.iterations
44
>>> solution.function_calls
46
``` You can also modify the solution accuracy. In this case we are willing to use a lower tolerance (up to the 2nd decimal place), to speed up the solution time: ```
>>> root, solution = bisect(jacketed_CSTR, 300, 310, args=(concA, Tin, Tj), full_output=True, xtol=1E-2)
>>> solution.converged
True
>>> solution.root
303.427734375
>>> solution.iterations
10
>>> solution.function_calls
12
``` All these methods find the root, \(T=303.4\)K. |

To use the Newton-Raphson method, you must code a second function that returns the derivative of \(f(x)\), i.e. \(f'(x)\). Then you can call it from Python as follows:

`root = newton(f, x_initial, f_derivative)`

where `f` and `f_derivative` are the names of the two functions and `x_initial` is your initial guess.