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: full_output=True >>> 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.