28.5: Advanced topics
-
- Last updated
- Save as PDF
This section introduces a few more advanced topics that allow you to use computer programming to simplifying many tasks. In this section, we will show you how you can write your own program to numerically estimate the value of an integral of any function.
Defining your own functions
Although Python provides many modules and functions, it is often useful to be able to define your own functions. For example, suppose that you would like to define a function that calculates \(\frac{1}{2}x^{2} +\frac{1}{4}x^{3} + \cos(2x)\), for a given value of \(x\). This is done easily using the
def
keyword in Python:
python code \(\PageIndex{1}\)
Defining a function
#import the math module in order to use cos import math as m #define our function and call it myfunction: def myfunction(x): return x**2 / 3 + x**3 / 4 + m.cos(2*x) #Test our function by printing out the result of evaluating it as x = 3 print ( myfunction(3) )
Output
10.710170286650365
A few things to note about the code above:
-
Functions are defined using the
def
keyword followed by the name that we choose for the function (in our case,myfunction
) -
If functions take arguments, those are specified in parenthesis after the name of the function (in our case, we have one argument that we chose to call
x
) - After the name of the function and the arguments, we place a colon
- The code that belongs to the function, after the colon, must be indented (this allows Python to know where the code for the function ends)
-
The function can “return” a value; this is done by using the
return
keyword. -
We used the “operator”
**
to take the power of a number (x**2
), and the operator*
, to multiply numbers. Python would not understand something like2x
; you need to use the multiplication operator, i.e.2*x
.
In the example above, we wrote a Python function to represent a mathematical function. However, one can write a function to execute any set of tasks, not just to apply a mathematical function. Python functions are very useful in order to avoid having to repeatedly type the same code.
Recall that the
numpy
module allows us to apply functions to arrays of numbers, instead of a single number. We can modify the code above slightly so that, if the argument to the function,
x
, is an array, the function will gracefully return an array of numbers to which the function has been applied. This is done by simply replacing the call to the
math
version of the
cos
function by using the
numpy
version:
python code \(\PageIndex{2}\)
Defining a function that works on an array
#import the numpy module in order to use cos to an array import numpy as np #define our function and call it myfunction: def myfunction(x): return x**2 / 3 +x**3 / 4 + np.cos(2*x) #Test our function by printing out the result and evaluating it at x = 3 (same as before) print( myfunction(3) ) #Test it with an array xvals = np.array([1, 2, 3]) print ( myfunction(xvals) )
Output
10.710170286650365 [ 0.1671865 2.67968971 10.71017029]
where we created the array
xvals
using the
numpy
module.
Using a loop to calculate an integral
The ability to define our own functions in Python allows us to easily simplify complex tasks. Using “loops” is another way that computer programming can greatly simplify calculations that would otherwise be very tedious. In a loop, one is able to repeat the same task many times. The example below simply prints out a statement five times:
Python code \(\PageIndex{3}\)
A simple loop
#A loop to print out a statement 5 times: for i in range(5): print("The value of i is ", i)
Output
The value of i is 0 The value of i is 1 The value of i is 2 The value of i is 3 The value of i is 4
A few notes on the code above:
-
The loop is defined by using the keywords
for ... in
-
The value after the keyword
for
is the “iterator” variable and will have a different value each time that the code inside of the loop is run (in our case, we called the variablei
) -
The value after the keyword
in
is an array of values that the iterator will take -
The
range(N)
function returns an array ofN
integer values between0
andN-1
(in our case, this returns the five values \(0,1,2,3,4\)) - The code to be executed at each “iteration” of the loop is preceded by a colon and indented (in the same way as the code for a function also follows a colon and is indented)
We now have all of the tools to evaluate an integral numerically. Recall that the integral of the function \(f(x)\) between \(x_{a}\) and \(x_{b}\) is simply a sum:
\[\begin{aligned} \int_{x_{a}}^{x_{b}}f(x)dx &= \lim_{\Delta x \to 0}\sum _{i=0}^{i=N-1}f(x_{i})\Delta x \\[4pt] \Delta x&=\frac{x_{b}-x_{a}}{N} \\[4pt] x_{i}&=x_{a}+i\Delta x \end{aligned}\]
The limit of \(∆x → 0\) is equivalent to the limit \(N → ∞\). Our strategy for evaluating the integral is:
- Define a Python function for \(f(x)\).
-
Create an array,
xvals
, of \(N\) values of \(x\) between \(x_{a}\) and \(x_{b}\). -
Evaluate the function for all those values and store those into an array,
fvals
. -
Loop over all of the values in the array
fvals
, multiply them by \(∆x\), and sum them together.
Let’s use Python to evaluate the integral of the function \(f(x) = 4x^{3}+3x^{2}+5\) between \(x = 1\) and \(x = 5\):
Python Code \(\PageIndex{4}\)
Numerical integration of a function
#import numpy to work with arrays: import numpy as np #define our function def f(x): return 4*x**3 + 3*x**2 + 5 #Make N and the range of integration variables: N - 1000 xmin = 1 xmax = 5 #create the array of values of x between xmin and xmax xvals = np.linspace(xmin, xmax, N) #evaluate the function at all those values of x fvals = f(xvals) #calculate delta x deltax = (xmax - xmin) / N #initialize the sum to be zero: sum = 0 #loop over the values fvals and add them to the sum for fi in fvals: sum = sum + fi*deltax #print the result: print("The integral between {} and {} using {} steps is {:.2f}".format(xmin, xmax, N, sum))
Output
The integral between 1 and 5 using 1000 steps is 768.42
One can easily integrate the above function analytically and obtain the exact result of \(768\). The numerical answer will approach the exact answer as we make \(N\) bigger. Of course, the power of numerical integration is to use it when the function cannot be integrated analytically.
Exercise \(\PageIndex{1}\)
What value of \(N\) should you use above in order to get within \(0.01\) of the exact analytic answer?
- Answer
-