Computing gradients

Solving Problems Numerically

for i in range(xlen):

for j in range(ylen):

f[i, j] = func(grid[0][i], grid[1][0,j])

plt.figure(figsize=(5, 5))

c = plt.contour(grid[0].flatten(), grid[1].flatten(), f)

# Compute and plot gradient of function grad = np.gradient(f, dx, dy)

plt.quiver(grid[0].flatten(), grid[1].flatten(), grad[0], grad[1])

plt.savefig('contour.png')

plt.close()

A contour plot of the function is shown below, with overlaid vectors representing the gradient:

[ 216 ]

Chapter 8

What just happened?

Computing the gradient is actually the easiest aspect of this example. Since we are using NumPy, we first need to define a two-dimensional grid of points. We use the ogrid function from NumPy, which does exactly that. We then used two nested for loops to compute the value of the function at each point in the grid. We used the contour function from Pyplot to display the level contours of the function, and used the gradient function to compute the gradient at each point of the grid. gradient is very similar to diff, except that it expects an N-dimensional array as its first argument. The next N arguments specify the grid spacing for each dimension of the array. Finally, we used the quiver plotting function from Pyplot to display a vector that represents the gradient at each point of the grid. We can see that the gradient is normal to the contour lines at each point.

Numericalintegration

Numerical integration (known in older literature as "quadrature") is another fundamental operation in numerical mathematics.

Timeforaction–numericalintegration

Let's start by using Sage functions to numerically integrate a symbolic function of one variable:

var('x')

f(x) = e^x * cos(x) f.show()

a = 0 b = 8

p = plot(f, (x, a, b)) p.show(figsize=(4, 3))

print("Integral of f(x) from {0} to {1}:".format(a,b)) print(" Analytical definite integral: {0}"

.format(f.integral(x, a, b).n()))

integral_value, tolerance, num_evals, error_code = \ f.nintegral(x, a, b)

print(" Using nintegral: {0}".format(integral_value))

# also nintegrate

integral_value, tolerance = numerical_integral(f, a, b)

print(" Using numerical_integral: {0}".format(integral_value))

[ 217 ]

Solving Problems Numerically

A plot of the function, along with the integration results, is shown below:

What just happened?

We defined a symbolic function f(x) and defined an interval on which to compute the definite integral. As a point of reference, we used the symbolic integration method integral

with limits to compute the definite integral symbolically. We used the method n (which is a shortcut for numerical_approx) to convert the symbolic result into a decimal number for comparison to the results of numerical integration.

Sage has two types of numerical integration. We used the numerical integration method nintegral (also called nintegrate) to perform the integration numerically. These methods call Maxima to perform the computation, and return a tuple with four elements:

[ 218 ]

Chapter 8

We repeated the calculation with numerical_integral, which uses the GNU Scientific Library. This function returns a tuple that contains the integral value and an estimate of the error. numerical_integral can be used to integrate functions that are defined as callable symbolic expressions or as Python functions. Both methods of numerical integration use algorithms that automatically adapt to the function that is being integrated. We can see that all three types of integration obtain equivalent results.

NumericalintegrationwithNumPy

Sometimes, we may need to integrate a function that cannot be defined by a Python function or a symbolic expression. This situation may occur if the function values come from experimental measurements or simulation results. For example, if we have an array of flux values over time, we can integrate the flux over time to obtain the total amount of substance. If the values of a function are stored in an array, we can use NumPy to perform the integration.

Timeforaction–numericalintegrationwithNumPy

We'll repeat the calculation from the previous example with NumPy. Execute the following code and see what happens:

import numpy as np def func(x):

return np.exp(x) * np.cos(x)

a=0; b=8 dx = 0.1

for i in range(4):

x = np.arange(a, b + dx, dx) f = func(x)

print("dx = {0} integral = {1}".format(dx, np.trapz(f,x))) dx = dx / 10

The output is shown in the following screenshot:

[ 219 ]

Chapter 8

freq = np.fft.fftfreq(t.shape[-1], d=dt)

spectrum = np.fft.fftshift(spectrum) phase = np.fft.fftshift(phase)

freq = np.fft.fftshift(freq)

plt.figure(figsize=(6,3)) plt.plot(freq, spectrum) plt.title('Magnitude')

plt.axis([-1.5, 1.5, spectrum.min(),spectrum.max()]) plt.savefig('spectrum.png')

plt.figure(figsize=(6,3)) plt.plot(freq, phase) plt.title('Phase')

plt.axis([-1.5,1.5, phase.min(),phase.max()]) plt.savefig('phase.png')

plt.close()

The following plots show the function, its magnitude spectrum, and its phase spectrum:

[ 221 ]

Solving Problems Numerically

What just happened?

We first generated a time-domain signal to work with. The arange function was used to create a series of sample times, and the sinc function was used to compute the sinc

function at those times. The first plot shows the time-domain function that we are going to analyse. We then used the fft function to compute the Fourier transform, which returns a NumPy array of complex values. It is customary to plot the frequency spectrum as magnitude and phase, so we used the absolute function to get the magnitude and the angle function to get the phase. The fftfreq function accepts an array of sample times and the time step, and uses them to compute the frequency values for plotting the frequency spectrum. The fftshift function was used to shift the magnitude and phase arrays and the frequency array so that the zero-frequency point was in the centre of the array. This is helpful for making plots.

[ 222 ]

Chapter 8

Windowfunctions

The sinc function that we used in the previous example is actually an approximation to the true sinc function, which extends from negative infinity to positive infinity. Because we used finite data to compute the Fourier transform, we introduced a distortion into the frequency spectrum. One way of reducing this distortion is by applying a window function to the signal before performing the transform. All the common window functions are built into NumPy.

Timeforaction–plottingwindowfunctions

Let's plot the window functions available in NumPy:

import numpy as np

import matplotlib.pyplot as plt plt.figure(figsize=(6,4)) plt.plot(np.bartlett(51), label="Bartlett") plt.plot(np.blackman(51), label="Blackman") plt.plot(np.hamming(51), label="Hamming") plt.plot(np.hanning(51), label="Hanning") plt.plot(np.kaiser(51,3), label="Kaiser") plt.legend(loc='best') plt.savefig('window_functions.png') plt.close()

The output is shown below:

[ 223 ]

Solving Problems Numerically

What just happened?

Calling the function hamming(51) returns an array containing 51 samples of the Hamming window function. The sample at the centre of the array corresponds to the centre of the window, where it reaches the value of 1. The other window functions work the same way, except that the Kaiser window also requires a shape parameter.

Haveagohero–usingwindowfunctions

Apply some of the window functions to the sinc function in the example and compute the

DFT to see how the window functions affect the frequency spectrum.

HINT: Generate an array of samples of the window function that has the same number of samples as the array containing the sinc function, and multiply the two arrays before performing the DFT.

Solvingordinarydifferentialequations

In the previous chapter, we looked at some tools for finding exact solutions to ordinary differential equations (ODEs). Sometimes, finding an exact solution is impossible or impractical. Sage also has a powerful set of tools for solving ordinary differential equations numerically.

Timeforaction–solvingafirst-orderODE

Let's start by solving a single, first-order, ordinary differential equation. We'll compare the exact solution to the numerical solution:

var('x, y')

y = function('y', x)

ode = diff(y, x) + 1 / x * y == x * y^2 sol = desolve(ode, y, ics=[10, 1]) sol.show()

exact_plot = plot(sol, (x, 0.1, 10), color='red', marker='.')

rk4_plot = desolve_rk4(ode, y, ics=[10, 1], output='slope_field', end_points=[0, 10])

show(exact_plot + rk4_plot, figsize=(4, 3))

The symbolic solution is shown in the following screenshot, and its graph is plotted on the same axes as the numerical solution and the slope field:

[ 224 ]

Chapter 8

What just happened?

As we explained in the previous chapter, we defined a symbolic ordinary differential equation and found the solution using desolve. The warning messages occur during the calculation of the slope field, and can be safely ignored. We then used the desolve_rk4 function to compute the solution numerically. As its name implies, the functions use the 4th order Runge-Kutta method to integrate a single first-order ordinary differential equation numerically (the function rk from Maxima's dynamics package is used to perform the integration). desolve_rk4 accepts the following keyword arguments:

[ 225 ]

Solving Problems Numerically

SolvingasystemofODEs

Sage also has tools to solve systems of ordinary differential equations numerically. We can use this capability to solve any higher-order ODE, since any higher-order ODE can be broken down into a system of first-order differential equations. For the next two examples, we will look at the van der Pol oscillator:

It can be written in terms of two first-order ODEs:

Timeforaction–solvingahigher-orderODE

Let's see if Sage can find a symbolic solution:

var('t')

x = function('x', t) y = function('y', t) u = 1.0

de1 = diff(x,t) - y == 0

de2 = diff(y,t) + x - u * y * (1 - x^2) == 0

Van_der_Pol = [de1, de2] desolve_system(Van_der_Pol, [x, y], ics=[0, 2, 0])

[ 226 ]

Chapter 8

If you run this code, Sage will return an error. It is unable to find a symbolic solution. Now, let's try to solve it numerically:

var('t, x, y') u = 1.0

Van_der_Pol = [y, -x + u * y * (1 - x^2)]

sol = desolve_system_rk4(Van_der_Pol, [x, y], ivar=t, ics=[0, 2, 0], end_points=[0, 20])

t = [i for i, j, k in sol] x_sol = [j for i, j, k in sol] y_sol = [k for i, j, k in sol]

# Plot results

import matplotlib.pyplot as plt plt.figure(figsize=(4, 3)) plt.plot(t, x_sol) plt.xlabel('t') plt.ylabel('y(t)') plt.subplots_adjust(bottom=0.15) plt.savefig('Van_der_Pol_rk4.png') plt.close()

# Limit cycle in the phase plane plt.figure(figsize=(4, 4)) plt.plot(x_sol, y_sol) plt.axis('scaled') plt.xlabel('x')

plt.ylabel('y') plt.savefig('Van_der_Pol_rk4_phase.png') plt.close()

[ 227 ]

Solving Problems Numerically

The results should look like this:

[ 228 ]

Solving Problems Numerically

num_points=1000)

interpolator = T.interpolate_solution()

plot(interpolator, (0, 20), axes_labels=('t','y'), figsize=(4,3))

The output is shown in the following screenshot:

What just happened?

This solver has a very different interface from the one in the previous example. The ordinary differential equations are defined by a Python function called f_1. The first argument to this function is the independent variable, the second is a list of dependent variables, and the third is a list of parameters. y[0] represents the variable y, and y[1] represents x. In this case, the parameter list has only one element, which is the value of u. We used similar syntax

to define the Jacobian function, j_1. The Jacobian is only required for certain integration algorithms; see the ode_solver documentation for more details.

The code T = ode_solver() creates a solver object. We used the attribute algorithm to choose the rk8pd (Runge-Kutta Prince-Dormand 8,9) algorithm, and used the attributes function and jacobian to pass in the functions we just defined. We called the method ode_solve to solve the system, which accepts the following arguments:

[ 230 ]

Solving Problems Numerically

c4_plot = implicit_plot(y == 0, (x,0,2), (y,0,2)) min_plot = point(sol['x'], color='red', size=50)

rp = region_plot([2 * x + y <= 3, x + 2 * y <= 3, x >= 0, y >= 0], (x,0,2), (y,0,2))

g = graphics_array([c1_plot+c2_plot+c3_plot+c4_plot+min_plot, rp], 1, 2)

g.show(aspect_ratio=1)

The output is shown below:

What just happened?

A linear program consists of a linear function to be minimized or maximized, and a set of linear constraints. In this case, there are only two variables, which makes it easier to visualize what's going on. We specified the function as a vector of coefficients called c. We created a matrix G that represents the left-hand side of the constraints. Each row in the matrix is a vector of coefficients for that constraint. Note that linear_program expects all constraints to specified as lhs < rhs, so we have to enter the third and fourth constraints with negative coefficients. The right-hand sides of the constraints were used to form another vector called h. The two vectors and the matrix were passed as arguments to linear_program, which returns the solution as a dictionary. The dictionary contains information about the solution, the minimized parameters, the slack variables, and the solution to the dual program (for more about slack variables and the dual program, see http://en.wikipedia.org/ wiki/Linear_programming or an optimization textbook). linear_program accepts

an optional argument with the keyword solver, which determines which solver is used. If the value is None (the default), the solver from the Python package CVXOPT is used. If solver='glpk', the solver from the GNU Linear Programming Kit (GLPK) is used.

[ 232 ]

Chapter 8

The plotting features of Sage can help us visualize the problem and the solution. We formulated the constraints as equalities, and utilized the implicit_plot function to plot the resulting lines in the x-y plane. We used the point function to add a dot to indicate the location of the minimum. We also used region_plot to plot the region defined by the inequalities. Going further, we could use a contour plot to show how the objective function varies over the domain.

Fittingafunctiontoanoisydataset

Fitting a function to data is one of the most common applications of numerical optimization for engineers and scientists. This is often referred to as "least squares fitting," because an optimal fit is achieved when the sum of squared errors is minimized.

Timeforaction–leastsquaresfitting

The following is adapted from an example in the Sage documentation:

var('a, b, c, x') set_random_seed(0.0)

data = [(i, 1.2 * sin(0.5 * i - 0.2) + 0.1 * normalvariate(0, 1)) for i in xsrange(0, 4 * pi, 0.2)]

data_plot = list_plot(data)

model(x) = a * sin(b * x - c)

fitted_params = find_fit(data, model, solution_dict=True) print("a = {0}".format(fitted_params[a]))

print("b = {0}".format(fitted_params[b])) print("c = {0}".format(fitted_params[c]))

g(x) = model.subs(a=fitted_params[a], b=fitted_params[b], c=fitted_ params[c])

fitted_plot = plot(g(x), (x, 0, 4 * pi), color='red') show(data_plot + fitted_plot, figsize=(4, 3))

[ 233 ]

Solving Problems Numerically

The noisy data, and the fitted function, are shown below:

What just happened?

After declaring some symbolic variables, we called the set_random_seed function to fix the seed for the pseudo-random number generator. A pseudo-random number generator will always generate the same sequence of numbers for a given seed value. Setting the seed ensures that the plot you create will be identical to the one in the example. We used a fairly complex list comprehension to construct the pseudo-random data set. We also defined a callable symbolic expression called model that has one variable and three variable parameters. We then used find_fit to vary the parameters until an optimal fit was achieved. find_fit has only two required arguments: a list containing the data and the model to be fitted. The keyword argument solution_dict=True tells the function to return the fitted values as a dictionary, rather than a list of relations. Having the parameter values in a dictionary allows us to substitute the values back into the model, using the subs method. Finally, we plotted the fitted function. find_fit accepts several additional optional arguments:

[ 234 ]