Only the first nine parts of the Basic Tutorial are included in this draft.

Basic Tutorial (part 2 of 10)

2 Numerical Calculations

In this second part of the Basic Tutorial you will learn about:

• Solving equations
• Finding minima and maxima
• Numerical calculus
• Numerical linear algebra
• Numbers
• Numerical error

Solving equations

The two most general functions in Mathematica for solving equations numerically are FindRoot and NSolve.  A third function, , is available for computing the roots of a single polynomial equation.  Systems of linear equations are an important special case. Functions for working with systems of linear equations are discussed later in this tutorial in the section on linear algebra.

The FindRoot function computes solutions for any equation or system of equations. Here is an example showing the use of FindRoot to find a root (a zero) of the function .  The first argument is the equation and the second argument is the list which gives the variable and a starting value for the variable.  The result is a list containing a rule that gives the value of the variable at the root.

Starting values for the variables in FindRoot are typically chosen by using estimates of the expected solutions, by trial and error, or by making plots. Here is a plot of the function showing the solution near 1.51213 that was returned by FindRoot.

The plot reveals a second solution when x is near 0.6. This solution can be computed by giving a different starting value for the variable.

In general, the choice of starting values in FindRoot has a strong effect both on the efficiency of the calculation and in determining whether or not FindRoot converges to the solution that you want.

FindRoot can also be used with two starting values for each variable.  With one starting value for each variable, FindRoot uses algorithms that require derivatives of the functions in the equations.  With two starting values for each variable, FindRoot uses a different algorithm (a secant method) for which derivatives are not required. Specifying two starting values for each variable is useful if calculation of the derivatives (the Jacobian of the equations) is difficult or requires a significant amount of time or memory.

Here is an example showing the use of FindRoot with two starting values for the variable.

FindRoot can be used with any number of equations. Here is an example showing the use of FindRoot to compute a solution for a pair of equations.  Starting values must be specified for each variable.  The solution is returned as a list of rules.

The NSolve function computes all of the solutions of a polynomial equation or system of polynomial equations, and sets of solutions for non-polynomial equations.  Here is an example showing the use of NSolve to find all four solutions of the pair of polynomial equations used in the preceding example.  The third solution in the result from NSolve corresponds to the solution returned by FindRoot.  The other solutions can be computed by giving different starting values for the variables in FindRoot.

You can get an equivalent result by applying N to the output of Solve.  In most cases, any difference between the result from NSolve and the result of applying N to the output of Solve will be small and are a consequence of numerical error.

Finding minima and maxima

The FindMinimum function is used for finding minima of a nonlinear function of one or more variables.  As with FindRoot, it is necessary to give starting values for the variables in FindMinimum.  

Here is an example showing the use of FindMinimum to find a minimum of the cubic polynomial .  As in FindRoot, the first argument in FindMinimum gives the function and the second argument is a list giving the variable and a starting value for the variable.  The result is a list giving the value of the function at the minimum and a rule that gives the value of the variable at the minimum.

A plot of this function shows the local minimum returned by FindMinimum and a nearby local maximum.

The local maximum of this function can be found by using FindMinimum to compute the local minimum of the negative of the function. With this approach, FindMinimum returns the negative of the value of the original function at the maximum and a rule giving the value of x at that point.

As might be expected from the plot, starting values must be chosen carefully in this example so that the search for a minimum does not enter a region in which the function decreases without bound. If that happens, FindMinimum will generate a warning message and the search for a minimum will fail.  For example:

The solution to this problem is to choose a better starting value for the variable.

The FindMinimum function can also be used with two starting values for each variable. With one starting value for each variable, FindMinimum uses algorithms that require derivatives (the gradient) of the first argument.  With two starting values for each variable, these derivatives are not needed.

Here is an example in which two starting values are specified for each variable. Note that the starting values do not need to bracket the solution. This example is sufficiently simple that almost any reasonable choice of starting values will work, including giving one starting value for each variable.

Numerical integration

Integrals can be computed numerically using the NIntegrate function. Here is an example showing the use of NIntegrate to compute the integral of Exp[x] for x from 0 to 1. The first argument in this example gives the integrand and the second argument gives the variable and the limits of integration.

The NIntegrate function is designed to handle integrable singularities at the endpoints of the range of integration automatically, as well as infinite ranges of integration. Here is an example in which the integrand is singular at the lower endpoint, and the range of integration is infinite.

Multiple integrals are specified by including additional arguments specifying the variables and the range of integration for each variable. The range specification for each variable can depend on variables that occur in earlier arguments. For example:

Nearly all difficulties with numerical integration can be resolved by considering the nature of the integrand and the limitations of numerical integration. Here, for example, is a calculation that returns a poor result because of the singularity in the range of integration. Although NIntegrate can handle singularities at the endpoints of the range of integration, all conventional numerical integration algorithms effectively work by approximating the integrand with smooth functions (usually polynomials), and will have difficulty if the integrand includes sharp features that cannot be approximated by smooth functions.

This problem can be resolved by providing the integration algorithm with information about the location of sharp features of the integrand. In NIntegrate, this can be done by including the location of the singularity as an additional limit of integration. With this change, the calculation is completed without difficulty.

Solving differential equations

Numerical solutions of differential equations can be computed in Mathematica using the NDSolve function. This function provides access to algorithms for solving initial value problems, boundary value problems, and partial differential equations.

An initial value problem is a calculation specified by one or more ordinary differential equations together with values for the solution and any necessary derivatives at a single value of the independent variable. The values for the solution are referred to generically as initial conditions.

Here is an example showing the use of NDSolve to compute the solution to the differential equation with initial condition . The notation is the short notation for Derivative[1][f][x], which denotes the first derivative of with respect to . The solution is assigned as the value of sol for use in subsequent calculations.

The first argument in NDSolve is a list that gives the differential equation and the initial condition. The second argument is the function for which to solve. The third argument is a list that gives the independent variable and the starting and ending values for that variable.

The result is expressed using an InterpolatingFunction expression. In general, the result from NDSolve will be a list of solutions, each of which is a list of rules, just as in results from Solve. In this example there is only one solution, and that solution consists of a single rule. The right side of that rule is an InterpolatingFunction expression.

The InterpolatingFunction expression can be used much like any other function. For example, if the InterpolatingFunction expression returned by NDSolve is substituted for f in f[5], the result will be the solution to the differential equation when x is 5. This substitution can be made by evaluating the following input:

This result is a list containing a number because the value of sol is a list that contains the solution. In this example there is only one solution, which can be selected using sol[[1]]. With this change the result of the substitution will be a number rather than a list containing a number.

The same substitution can be used to make a plot of the solution.

Here is a more complicated initial value problem involving a system of three ordinary differential equations. This example also has two solutions, so the result is a list with two elements, each of which is a list of three rules.

An boundary value problem differs from an initial value problem in that values for the solution in a boundary value problem are specified at more than one value of the independent variable. Values for the solution in a boundary value problem are called boundary conditions. Here is a boundary value problem in which the value of the solution f[x] is specified by the boundary conditions and . The differential equation in this example could describe the height of a projectile, in which case the boundary conditions might indicate that the projectile should be on the ground when x is 0 and when x is 1.

NDSolve can also be used to solve partial differential equations. A partial differential equation is a differential equation with more than one independent variable. In Version 4 of Mathematica, the numerical partial differential equation solver in NDSolve can handle partial differential equations with two independent variables.

Here is an example showing the use of NDSolve to solve a problem involving a common partial differential equation called the wave equation. Derivatives are entered here using FullForm notation. Notations such as f'[x] that can be used for entering ordinary differential equations are unambiguous only for specifying derivatives of functions of a single variable.

The equations in this example describe the oscillations of a string of length Pi, where gives the displacement of the string at point x and time t. The equation gives the initial displacement (the value of when t is 0), and the equations and specify that the string is fixed at the endpoints (the value of is zero for all values of t when x is 0 and when x is Pi).

The equation is the wave equation. denotes the second derivative with respect to the first argument of f evaluated when the first argument is x and the second argument is t. denotes the second derivative with respect to the second argument of f evaluated when the first argument is x and the second argument is t. The equation specifies that the first derivative of f with respect to the second variable is zero for all x when t is zero (the string is initially motionless). The meaning of these notations involving Derivative is discussion in the documentation for Derivative.

The result from NDSolve is again given in terms of an InterpolatingFunction expression. A plot of this result can be constructed using the Plot3D function. The Plot3D function is discussed in Part 4 of this tutorial.

Evaluation in numerical functions

Many numerical functions hold their arguments unevaluated until numerical values are assigned to the variable or variables of the function. Of the functions that have been discussed here so far, FindRoot, FindMinimum, and NIntegrate all have this property. Hold arguments unevaluated enables these functions to work even when the variables have been assigned values outside of the function, and provides for purely numerical calculations without the symbolic calculation that would otherwise occur before the function is called.

In most cases this detail of the design of these functions can be ignored. There are, however, two common situations in which important to be aware of how these functions evaluate their arguments.

One situation occurs when evaluation of the function arguments changes the structure of the arguments. For example, the following input to NIntegrate fails because the symbol lims in the second argument is held unevaluated, and so is not recognized as providing valid limits of integration.

There are several ways to address this type of problem, perhaps the most common of which is to use Evaluate to force evaluation of the function arguments.

A second situation in which evaluation of function arguments is a concern occurs for calculations that can, at least in principle, be done using a mixture of symbolic operations and numerical operations.

To illustrate this effect, consider the following three methods of combining Integrate and NIntegrate to evaluate a double integral, where in each method Integrate is used to evaluate the inner integral symbolically and NIntegrate is used to evaluate the outer integral numerically.

In the first calculation, since NIntegrate holds the integrand unevaluated and evaluates the integrand for each new numerical value that it assigns to the variable, the symbolic integral in the integrand is evaluated repeatedly for each new value of the variable. This is a relatively inefficient approach.

In the second calculation, Evaluate is used to force evaluation of the symbolic integral for a symbolic value of x, so that the symbolic integral is only evaluated once.

In the third calculation the symbolic integral is evaluated outside of NIntegrate and assigned as the value of expr, which is then used as the integrand in NIntegrate. With this approach the symbolic integral is again only evaluated once.

Although efficiency is improved in this example by forcing a symbolic calculation, this is not the case in all examples. In many examples, it is more efficient to do several numerical calculations rather than a single symbolic calculation. The best strategy is to analyze each calculation and choose the optimal approach for that calculation.

The feature of holding arguments unevaluated in functions such as NIntegrate is controlled by the HoldAll attribute. Attributes are properties of symbols that control various behaviors, such as whether or arguments are held unevaluated when a symbol is used as a function. You can use the Attributes function to check the attributes of a symbol. The HoldAll attribute is include in the attributes of NIntegrate.

For more information on attributes, see section 2.5.3 in The Mathematica Book.

Other numerical operations on functions

There are many other numerical operations on functions that are included with Mathematica.  For more information see Section 3.9 of the Mathematica book or the section on Numerical Mathematics Packages in the documentation for standard add-on packages.

Vectors, matrices, and tensors

For a basic discussion of vectors and matrices, see Section 1.8.3 in The Mathematica Book. Vectors, matrices, and tensors are represented in Mathematica using lists. A list of elements is a vector, and a list of lists is a matrix. More deeply nested lists represent other tensors.

Here is a vector with three components, assigned as the value of v.

A matrix is a list in which all of the elements are lists of the same length. Here is a matrix, assigned as the value of m.

Multiplication of a vector or a matrix by a constant can be done using ordinary multiplication. For example:

Vector addition is done using ordinary addition. For example:

Ordinary addition of a constant and a vector or a constant and a matrix adds the constant to each element of the vector or matrix. For example:

Inner products (dot products) of two vectors can be computed using the Dot function. Dot products can be entered using a period (.) between the operands. For example, Dot[v,v] is equivalent to .

Dot is also used to compute the dot product of a vector and a matrix, or the dot product of a matrix and a matrix. For example:

Ordinary multiplication of a vector by another vector, or a matrix by another matrix, does term by term multiplication of the elements of the vector, or term by term multiplication of the elements in the matrix. This is not the same as a dot product. For example:

Numerical linear algebra

Operations with vectors and matrices are referred to generically as operations in linear algebra, in part because these operations are often used for solving systems of linear equations. For example, the set of three linear equations in the input:

can be represented using matrices and solved using LinearSolve:

The first argument in LinearSolve is a matrix that gives the coefficients of the variables on the left side of each equation, and the second argument is a vector that gives the constant on the right side of each equation.

Mathematica includes a broad assortment of functions for common matrix operations. For example, the Inverse function provides another way of solving this system of linear equations.

For more information about functions for working with vectors and matrices in Mathematica, see Section 3.7 in The Mathematica Book.

Nearly all of the matrix functions in Mathematica can in principle be used both for numerical matrices (matrices in which all of the matrix elements are numbers) or for symbolic matrices (matrices in which some or all of the matrix elements contain symbols). However, although operations with symbolic matrices may be conceptually similar to operations with numerical matrices, operations with symbolic matrices almost always use different algorithms that operations with numerical matrices, and usually require significantly more computer time and memory.

The difference can be dramatic. For example, although LinearSolve can compute a solution for 100 simultaneous linear equations with numerical coefficients in a fraction of a second, the analogous calculation with general symbolic matrices would vastly exceed the capacity of any realistic computer.

As a simple illustration of this difference in computational requirements between symbolic linear algebra and numerical linear algebra, note that the result for the inverse of a general symbolic matrix is much larger than the result for the inverse of a numerical matrix. This difference becomes much greater for larger matrices. For matrices larger than about , the symbolic result is so large that storing or manipulating that result will exhaust the memory of most computers.

Although symbolic linear algebra can of course be very useful, it is also useful to keep in mind that some symbolic calculations can be prohibitively large, and to be careful in setting up numerical calculations to avoid doing symbolic calculations by mistake.

Many operations with vectors and matrices are handled in Mathematica using the same functions that are used for manipulating other types of expressions. To demonstrate some of these operations, here is matrix assigned as the value of m, and displayed using MatrixForm.

This picks out an element in the matrix m.

This picks out the second row in m.

This picks out the second column in m.

This picks out the block in the lower right corner of m.

This assigns a new value to an element of m. The effect of this assignment can be seen most easily by looking at the new value of m.

This divides the second row of the matrix by the value of the first element in that row, subtracts that result from the third row, and assigns that difference as the new value of the third row.

There are many other operations that are useful both for manipulation of general expressions and for common operations in linear algebra.

Exact and inexact numbers

Numbers are of course a basic element in all numerical calculations. All numbers are either exact or inexact. An exact number is a number like or that is represented exactly, with no numerical uncertainty.  An inexact number is a number like that is specified by giving some finite number of digits.  The digits that fall beyond the specified digits in an inexact number are a source of numerical error.

In simple calculations, exact inputs lead to exact results, and inexact inputs lead to inexact results.  For example, if you evaluate , Mathematica will return the exact result , and if you evaluate , using the inexact number rather than the exact number , Mathematica will return an inexact approximation for .

Most computer programs, such as the programs found on a simple pocket calculator, are not designed to represent exact quantities like or .  If you ask a simple pocket calculator for the square root of 3, for example, you will normally get an inexact approximation that gives the first ten or twenty digits of the square root of 3.  If you evaluate in Mathematica the result will be the exact quantity .

Although exact numbers avoid a significant source of numerical error, most numerical calculations are for practical reasons done using inexact numbers.

You can get inexact results in Mathematica either by starting with inexact inputs or by using a function that is designed to give inexact results.  One function that gives inexact results is FindRoot, which was encountered in the first part of this tutorial.  The FindRoot function will give inexact results even if the input is exact. For example:

Another function that gives inexact results is the N function, which is often used to compute inexact approximations to exact expressions.  For example, you can use to compute an inexact approximation for

Other functions that give inexact results even if the inputs are exact include NSolve, FindMinimum, NSum, NProduct, NIntegrate, and NDSolve.

Precision and accuracy

In Mathematica, the precision of an inexact number is defined as the number of significant digits in the number, and the accuracy of an inexact number is defined as the position of the last significant digit relative to the decimal point.  For example, the number , with 27 digits, has a precision of 27.  The accuracy of this number is 24, since the last significant digit is 24 digits after the decimal point.  Precision and accuracy can also be negative.  If the last significant digit occurs to the left of the decimal point, for example, then the accuracy will be negative.

The Precision function can be used to get the precision of a number, and the Accuracy function can be used to get the accuracy of a number.  If is an expression, gives the precision of the least precise number in and gives the accuracy of the least accuracy number in .

Mathematica supports both machine numbers and variable-precision numbers.  A machine number is a number that is represented directly using a corresponding numeric data type on your computer.  Variable-precision numbers are stored and manipulated in software within Mathematica.  Arithmetic with machine numbers is usually handled by the hardware of your computer.

Most numerical calculations, including most default numerical calculations in Mathematica, are done using machine numbers.  Machine numbers use less memory that variable-precision numbers, and calculations with machine numbers are usually faster than calculations with variable-precision numbers.  Variable-precision numbers are used primarily for high-precision calculations and in calculations that require greater control over numerical error.

Machine numbers are fixed-precision numbers.  All machine numbers in a given version of Mathematica have the same number of digits.  The value of the variable $MachinePrecision gives the precision of machine numbers in your version of Mathematica.

The precision of a variable-precision number can be any real number between the values of $MaxPrecision and $MinPrecision.

Calculations in applied mathematics rarely require more than a few hundred digits of precision.  For precision greater than the value of $MaxPrecision you can set $MaxPrecision to a larger value.  Calculations with more than a few million digits of precision can require a great deal of time and memory, and are usually of interest only in special applications.

The precision of variable-precision results is adjusted automatically to include only the digits that are numerically meaningful.  Here is an example of a subtraction in which many of the digits cancel out, so that the precision of the result is considerably less than the precision of the inputs.

If the corresponding calculation were to be done using fixed-precision arithmetic, the number of digits in the result would be the same as the number of digits in the input, even though many of those digits would not necessarily be correct.  This issue is discussed in greater detail in the section below on numerical error.

Entering and displaying numbers

Numbers are usually entered in to Mathematica using decimal digits.  Inexact numbers are distinguished from exact numbers by the presence of a decimal point.

Numbers can be entered in scientific notation by using *^ to identify the exponent.  For example, 3.14159*^25 indicates multiplied by .

A number in input is assumed to be a machine number if the number of digits is less than or equal to the value of $MachinePrecision and if the magnitude of the number falls between the value of $MinMachineNumber and the value of $MaxMachineNumber.  The value of $MinMachineNumber gives the smallest positive number that can be represented as a machine number on your computer, and the value of $MaxMachineNumber gives the largest number that can be represented as a machine number on your computer.

Numbers outside of this range, and numbers with precision greater than the value of $MachinePrecision, are represented using variable-precision numbers.

You can use the function MachineNumberQ to find out if a number is an inexact machine number.   returns True if is an inexact machine number, and False otherwise.  For example:

In output, machine numbers are displayed by default using 6 digits, with trailing zeroes omitted.  Machine numbers that would require more than 6 digits either to the right or to the left of the decimal point are shown in scientific notation.

Here are some examples showing the default display of machine numbers.

Variable-precision numbers are displayed by default by giving all significant digits, including trailing zeroes.  Variable-precision numbers larger than or equal to , or smaller than or equal to , are shown in scientific notation.  For example:

A useful function for displaying numbers is InputForm, which uses a format that can be read back in to Mathematica to reproduce the original number.  The notations used for numbers in InputForm format are described in Section A.2.5 of The Mathematica Book.

InputForm notation can be used to show all of the non-zero digits in a machine number.

3.14159255358979

InputForm notation can also be used in input to specify the precision or accuracy of a variable-precision number.  For example:

3.14159`30

Other methods for controlling the display of numbers are described in Section 2.8.7 of The Mathematica Book, and in the notes for functions such as NumberForm, ScientificForm, PaddedForm, and FortranForm.

This section provides a brief introduction to the subject of numerical error.  Numerical error is a fundamental topic in numerical analysis. An appreciation for this topic is essential in nearly all numerical work.

Numerical error is inevitable in all calculations with inexact numbers.  Even in simple examples numerical error can dominate a result and can do so with no obvious indication that anything has gone wrong.  Fortunately it is relatively easy to learn to anticipate and allow for the consequences of numerical error.

Numerical errors can be classified as roundoff errors and truncation errors.  Most numerical calculations are affected by both types of error.

Roundoff error is error associated with using limited precision in calculations.  Suppose, for example, that a calculation is to be done using numbers with sixteen decimal digits of precision.  Since most numbers cannot be represented exactly using sixteen decimal digits, numerical inputs and results must be rounded to the best available sixteen-digit approximations.  For example, the best available sixteen-digit approximation for the number is , which is slightly smaller than .  The difference between and is a form of roundoff error. In a typical calculation, roundoff errors in inputs and in intermediate results will propagate through the calculation and will lead to a corresponding roundoff error in the final result.

Roundoff error can be reduced either by using higher precision or by rearranging the calculation to stabilize it against the effects of roundoff error. Use of high precision is a simple and straightforward approach to reducing roundoff error, but since high-precision calculations require more time and memory than low-precision calculations, this approach is not appropriate in all applications. In most Mathematica functions the precision used for computing intermediate results is controlled by the WorkingPrecision option. Rearranging the calculation to improve numerical stability is the preferred approach if computational efficiency is a dominant concern. This approach can be difficult to implement, and is beyond the scope of this tutorial. This topic is covered extensively in books and classes on numerical analysis. For a classic example of rearranging a calculation to improve numerical stability, see the discussion of the Horner form of a polynomial in the documentation for the standard Algebra`Horner` package.

The second type of numerical error is truncation error.  Truncation error is the error that results from doing a limited number of steps in a calculation that in principle requires an infinite (or very large) number of steps.  Numerical approximations for integrals, for example, are computed by evaluating the integrand at a finite number of points and using those points to estimate the value of the integral. The difference between this estimate and the actual value of the integral is a form of truncation error.  Truncation error affects many numerical calculations, including numerical calculation of roots, numerical minimization, numerical integration, and numerical solution of differential equations.

Truncation error can be reduced either by increasing the number of steps in the calculation or by choosing an algorithm that is better suited to a particular calculation. In most cases the number of steps is controlled indirectly, by setting a tolerance for the error in the result. The  number of steps is then increased automatically until an estimate of the error falls below the specified tolerance, or until some limiting number of steps is reached. In most Mathematica functions, tolerances are specified using options such as AccuracyGoal and PrecisionGoal, and limits on the number of steps are specified using options such MaxSteps, MaxRecursion, or MaxIterations. The other approach to reducing truncation error is to choose a different algorithm. Strategies for choosing algorithms are again beyond the scope of this tutorial, but are discussed in books and classes on numerical computation. The choice of algorithm in most Mathematica functions is controlled either by the form of the input or by the setting of the Method option for the function. In most examples Mathematica will choose an optimal (or nearly optimal) algorithm automatically.

Before taking steps to reduce numerical error it is important to do at least some simple analysis to determine the origin of the error. The is little reason to increase precision to reduce roundoff error, for example, if the primary source of error is truncation error. This topic will be revisited in the discussion in the next section on high-precision calculations.

It is also worth noting that almost any numerical calculation can give wildly incorrect results, and can do so without warning. For example, if an important feature of an integrand happens to fall between the points that are sampled by a numerical integration algorithm, the effect of that feature will be omitted in the result. It is not possible by purely numerical methods to avoid this type of error. For this reason it is useful to develop at least a passing understanding of the limitations of numerical algorithms before making extensive use of these functions.

High-precision calculations

There are two relatively obvious considerations to keep in mind for any high-precision calculation.

First, high-precision calculations are normally of interest only if all inputs and all intermediate calculations are done using high precision. There is usually little reason, for example, to do a high-precision calculation with low-precision inputs.

The second consideration is that using high precision to reduce roundoff error should in most cases be accompanied by adjusting the tolerances that control truncation error.

To illustrate these points, here is an example showing the use of FindMinimum to compute a local minimum of a quartic polynomial. By default, this calculation is done using machine arithmetic, and gives machine numbers in the result.

The FindMinimum function has seven options, most of which are related to controlling numerical error.

Increasing the values of PrecisionGoal and AccuracyGoal will not give a high-precision calculation because the default value of the WorkingPrecision option is such that intermediate results are computed using machine precision.

Setting the WorkingPrecision option to a higher value, however, does give a high-precision calculation, because the AccuracyGoal and PrecisionGoal options are automatically reset to comparable values.

If the input includes a low-precision number, FindMinimum will generate a warning message indicating that the result may not be correct. In this example, the precision of the machine number 3.0 in the input is too low to give a reliable calculation at the precision specified by the WorkingPrecision option. Only about the first sixteen digits in this result are correct.

Exercises

Exercise 1.1 Numerical approximations

Use FindRoot to compute the largest root of , and then compare that root with the exact result. You can get the exact result using Solve.

You can get a starting value to use in FindRoot either by plotting this polynomial or by using N to compute numerical approximations for the exact results from Solve. These numerical approximations can also be useful for identifying which of the three solutions returned by Solve corresponds to the largest root. The roots are not necessarily returned in numerical order.

Exercise 1.2 Starting values in FindMinimum

The function has several local minima. Use FindMinimum to find the local minimum with the largest value of x. The value of this function at that minimum is approximately 108.429.

Exercise 1.3 High-precision numerical integration

Here is an example that shows an attempt to get 50 digits of precision in the result of a numerical integral. This calculation fails because both the working precision and the precision of the integrand are too low.

Change this example to correctly give a result with 50 digits of precision.

Exercise 1.4 Plotting a solution computed by NDSolve

Make a plot of this solution from NDSolve over the full range of the solution.

Exercise 1.5 Fitting a curve to data

Here are two lists giving values for the independent and dependent (response) values in a data set. The value of xdata gives the values of the independent variable, and the value of ydata gives values for the response variable.

Use Fit to fit a linear combination of the functions {1,x,x^2} to these data. Transpose is useful for getting these data into the form required by Fit.

Exercise 1.6 Eigenvalues of a matrix

Here is a relatively inefficient way of computing the value of x for which the largest eigenvalue of the matrix m reaches a minimum.

Demonstrate a more efficient way of doing this calculation. The method shown above is inefficient because it involves constructing a large symbolic expression for the eigenvalues of this matrix. There are several simple ways of rearranging this calculation that avoid computing eigenvalues until after numerical values are inserted for x.

The semicolon at the end of the definition of expr is used to suppress the display of the expression for the eigenvalues. You can see that result by evaluating Eigenvalues[m].

The Max function in the definition of expr is used to pick out the largest eigenvalue. The Re function is included to discard tiny imaginary parts that arise as an artifact of numerical error, and that would interfere with identifying the maximum eigenvalue. This use of Re will probably be unnecessary in your solution.

You can use the Timing function to check the speed of your calculation. Since all of these calculations finish in a fraction of a second, you may be able to get more reliable timing information by using a Do loop to time many evaluations of each calculation. For example:

Exercise 1.7 Manipulating a matrix

With one use of the Part function, construct a matrix consisting of the elements in the corners of this matrix.

Exercise 1.8 Formatting numbers in a table

Consider the following example in which results from NDSolve are displayed in a table. The second column in this table is formatted using PaddedForm so that each entry in that column has four digits to the right of the decimal point.

0.`
0.5`
1.`
1.5`
2.`
2.5`
3.`
3.5`
4.`
4.5`
5.`

Change this table to so that the entries in the first column are formatted with one digit to the right of the decimal point. The resulting table should look like this:

Exercise 1.9 Resampling using interpolation

Many calculations require data sampled at regular intervals or on a regular grid. If the original data include missing entries or entries sampled at points with irregular spacing, interpolation or extrapolation can be used to generate data at regular points.

Here is a list of pairs representing data sampled at an irregular set of points.

Use Interpolation to resample these data at points from 0 to 100 in steps of 2. This can be done by evaluating Interpolation[data] to construct an InterpolatingFunction expression, and then using Table to construct a list of elements {k,f[k]} for k ranging from 0 to 100 in steps of 2, where f represents the InterpolatingFunction expression.

Plot the resampled data using ListPlot and compare that plot with the plot of the original data.

Exercise 1.10 Numerical differentiation

One of many useful functions in the standard NumericalMath add-on packages is ND, which is defined in the NumericalMath`NLimit` package. The ND function is used for computing derivatives numerically. Numerical differentiation is of interest most frequently for differentiating expressions that cannot be differentiated symbolically using the D function.

One common function that is not differentiated symbolically using D is the absolute value function Abs. The Abs function is not differentiated by D because the derivatives of Abs depend on direction. Mathematically, a function is said to be differentiable if the derivative exists and is the same in all directions. Since derivatives of Abs are different in different directions, Abs is not a differentiable function.

It is, however, possible to use ND to compute directional derivatives of expressions involving Abs. Direction is specified in ND using the Scale option. A positive value for the Scale option indicates a derivative from the right, and a negative value indicates a derivative from the left.

Use ND to plot the derivative of Exp[-Abs[x]] for x from -3 to 3. A separate plotting variable (say, x0) will be necessary in this calculation to distinguish the plotting variable from the variable of differentiation. The option Scale->Sign[x0] can be used so the derivative will be computed from the left for negative values of x0 and from the right for positive values of x0.