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.
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 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]
Similar to the MATLAB case above, what if we want to freeze the last 3 inputs, and create a new function that just accepts one input, the tank temperature?
from functools import partial
>>> concA = 0.4 # mol/m^3
>>> Tin = 300 # K
>>> Tj = 350 # K
>>> heat_balance = partial(jacketed_CSTR, CA=concA, Tin=Tin, Tj=Tj)
>>> heat_balance(T=305)
-0.077055755255798275
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 function available, called fzero. This 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.
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 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 *
# Create a partial function, as shown above:
>>> from functools import partial
>>> concA = 0.4 # mol/m^3
>>> Tin = 300 # K
>>> Tj = 350 # K
>>> heat_balance = partial(jacketed_CSTR, CA=concA, Tin=Tin, Tj=Tj)
# Apply the bisection method:
>>> bisect(heat_balance, 300, 310)
303.42075240691372
# Apply the Secant method: use the ``newton`` function, but don't give it f'(x)
>>> newton(heat_balance, 300)
# Apply Brent's method (similar to fzero in MATLAB):
>>> brentq(heat_balance, 300, 310)
303.42075240691321
|