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.

% 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]

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 use function handles in plot commands. Try this:

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

and use the plot to graphically locate the function's zero.

Python details: functions are objects

In Python, everything is an object. 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 a variable'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'>

Find zeros of a nonlinear function

• scipy.optimize.brentq
• scipy.optimize.newton
• brenth,
• ridder,
• bisect,
• newton
• fsolve
• fixed_point