Chapter 17 Objectives Recognizing that NewtonCotes integration formulas

Chapter 17 Objectives • Recognizing that Newton-Cotes integration formulas are based on the strategy of replacing a complicated function or tabulated data with a polynomial that is easy to integrate. • Knowing how to implement the following single application Newton-Cotes formulas: – Trapezoidal rule – Simpson’s 1/3 rule – Simpson’s 3/8 rule • Knowing how to implement the following composite Newton-Cotes formulas: – Trapezoidal rule – Simpson’s 3/8 rule

Objectives (cont) • Recognizing that even-segment-odd-point formulas like Simpson’s 1/3 rule achieve higher than expected accuracy. • Knowing how to use the trapezoidal rule to integrate unequally spaced data. • Understanding the difference between open and closed integration formulas.

Integration • Integration: is the total value, or summation, of f(x) dx over the range from a to b:

Newton-Cotes Formulas • The Newton-Cotes formulas are the most common numerical integration schemes. • Generally, they are based on replacing a complicated function or tabulated data with a polynomial that is easy to integrate: where fn(x) is an nth order interpolating polynomial.

Newton-Cotes Examples • The integrating function can be polynomials for any order - for example, (a) straight lines or (b) parabolas. • The integral can be approximated in one step or in a series of steps to improve accuracy.

The Trapezoidal Rule • The trapezoidal rule is the first of the Newton. Cotes closed integration formulas; it uses a straight-line approximation for the function:

Error of the Trapezoidal Rule • An estimate for the local truncation error of a single application of the trapezoidal rule is: where is somewhere between a and b. • This formula indicates that the error is dependent upon the curvature of the actual function as well as the distance between the points. • Error can thus be reduced by breaking the curve into parts.

Composite Trapezoidal Rule • Assuming n+1 data points are evenly spaced, there will be n intervals over which to integrate. • The total integral can be calculated by integrating each subinterval and then adding them together:

MATLAB Program

Simpson’s Rules • One drawback of the trapezoidal rule is that the error is related to the second derivative of the function. • More complicated approximation formulas can improve the accuracy for curves - these include using (a) 2 nd and (b) 3 rd order polynomials. • The formulas that result from taking the integrals under these polynomials are called Simpson’s rules.

Simpson’s 1/3 Rule • Simpson’s 1/3 rule corresponds to using second-order polynomials. Using the Lagrange form for a quadratic fit of three points: Integration over the three points simplifies to:

Error of Simpson’s 1/3 Rule • An estimate for the local truncation error of a single application of Simpson’s 1/3 rule is: where again is somewhere between a and b. • This formula indicates that the error is dependent upon the fourth-derivative of the actual function as well as the distance between the points. • Note that the error is dependent on the fifth power of the step size (rather than the third for the trapezoidal rule). • Error can thus be reduced by breaking the curve into parts.

Composite Simpson’s 1/3 Rule • Simpson’s 1/3 rule can be used on a set of subintervals in much the same way the trapezoidal rule was, except there must be an odd number of points. • Because of the heavy weighting of the internal points, the formula is a little more complicated than for the trapezoidal rule:

Simpson’s 3/8 Rule • Simpson’s 3/8 rule corresponds to using third-order polynomials to fit four points. Integration over the four points simplifies to: • Simpson’s 3/8 rule is generally used in concert with Simpson’s 1/3 rule when the number of segments is odd.

Higher-Order Formulas • Higher-order Newton-Cotes formulas may also be used - in general, the higher the order of the polynomial used, the higher the derivative of the function in the error estimate and the higher the power of the step size. • As in Simpson’s 1/3 and 3/8 rule, the even-segment -odd-point formulas have truncation errors that are the same order as formulas adding one more point. For this reason, the even-segment-odd-point formulas are usually the methods of preference.

Integration with Unequal Segments • Previous formulas were simplified based on equispaced data points - though this is not always the case. • The trapezoidal rule may be used with data containing unequal segments:

Integration Code for Unequal Segments

MATLAB Functions • MATLAB has built-in functions to evaluate integrals based on the trapezoidal rule • z = trapz(y) z = trapz(x, y) produces the integral of y with respect to x. If x is omitted, the program assumes h=1. • z = cumtrapz(y) z = cumtrapz(x, y) produces the cumulative integral of y with respect to x. If x is omitted, the program assumes h=1.

Multiple Integrals • Multiple integrals can be determined numerically by first integrating in one dimension, then a second, and so on for all dimensions of the problem.

Chapter 18 Objectives • Understanding how Richardson extrapolation provides a means to create a more accurate integral estimate by combining two less accurate estimates. • Understanding how Gauss quadrature provides superior integral estimates by picking optimal abscissas at which to evaluate the function. • Knowing how to use MATLAB’s built-in functions quad and quadl to integrate functions.

Richardson Extrapolation • Richard extrapolation methods use two estimates of an integral to compute a third, more accurate approximation. • If two O(h 2) estimates I(h 1) and I(h 2) are calculated for an integral using step sizes of h 1 and h 2, respectively, an improved O(h 4) estimate may be formed using: • For the special case where the interval is halved (h 2=h 1/2), this becomes:

Richardson Extrapolation (cont) • For the cases where there are two O(h 4) estimates and the interval is halved (hm=hl/2), an improved O(h 6) estimate may be formed using: • For the cases where there are two O(h 6) estimates and the interval is halved (hm=hl/2), an improved O(h 8) estimate may be formed using:

The Romberg Integration Algorithm • Note that the weighting factors for the Richardson extrapolation add up to 1 and that as accuracy increases, the approximation using the smaller step size is given greater weight. • In general, where ij+1, k-1 and ij, k-1 are the more and less accurate integrals, respectively, and ij, k is the new approximation. k is the level of integration and j is used to determine which approximation is more accurate.

Romberg Algorithm Iterations • The chart below shows the process by which lower level integrations are combined to produce more accurate estimates:

MATLAB Code for Romberg

Gauss Quadrature • Gauss quadrature describes a class of techniques for evaluating the area under a straight line by joining any two points on a curve rather than simply choosing the endpoints. • The key is to choose the line that balances the positive and negative errors.

Gauss-Legendre Formulas • The Gauss-Legendre formulas seem to optimize estimates to integrals for functions over intervals from -1 to 1. • Integrals over other intervals require a change in variables to set the limits from -1 to 1. • The integral estimates are of the form: where the ci and xi are calculated to ensure that the method exactly integrates up to (2 n-1)th order polynomials over the interval from -1 to 1.

Adaptive Quadrature • Methods such as Simpson’s 1/3 rule has a disadvantage in that it uses equally spaced points - if a function has regions of abrupt changes, small steps must be used over the entire domain to achieve a certain accuracy. • Adaptive quadrature methods for integrating functions automatically adjust the step size so that small steps are taken in regions of sharp variations and larger steps are taken where the function changes gradually.

Adaptive Quadrature in MATLAB • MATLAB has two built-in functions for implementing adaptive quadrature: – quad: uses adaptive Simpson quadrature; possibly more efficient for low accuracies or nonsmooth functions – quadl: uses Lobatto quadrature; possibly more efficient for high accuracies and smooth functions • q = quad(fun, a, b, tol, trace, p 1, p 2, …) – fun : function to be integrates – a, b: integration bounds – tol: desired absolute tolerance (default: 10 -6) – trace: flag to display details or not – p 1, p 2, …: extra parameters for fun – quadl has the same arguments

Chapter 19 Objectives • Understanding the application of high-accuracy numerical differentiation formulas for equispaced data. • Knowing how to evaluate derivatives for unequally spaced data. • Understanding how Richardson extrapolation is applied for numerical differentiation. • Recognizing the sensitivity of numerical differentiation to data error. • Knowing how to evaluate derivatives in MATLAB with the diff and gradient functions. • Knowing how to generate contour plots and vector fields with MATLAB.

Differentiation • The mathematical definition of a derivative begins with a difference approximation: and as x is allowed to approach zero, the difference becomes a derivative:

High-Accuracy Differentiation Formulas • Taylor series expansion can be used to generate high-accuracy formulas for derivatives by using linear algebra to combine the expansion around several points. • Three categories for the formula include forward finite-difference, backward finitedifference, and centered finite-difference.

Forward Finite-Difference

Backward Finite-Difference

Centered Finite-Difference

Richardson Extrapolation • As with integration, the Richardson extrapolation can be used to combine two lower-accuracy estimates of the derivative to produce a higher-accuracy estimate. • For the cases where there are two O(h 2) estimates and the interval is halved (h 2=h 1/2), an improved O(h 4) estimate may be formed using: • For the cases where there are two O(h 4) estimates and the interval is halved (h 2=h 1/2), an improved O(h 6) estimate may be formed using: • For the cases where there are two O(h 6) estimates and the interval is halved (h 2=h 1/2), an improved O(h 8) estimate may be formed using:

Unequally Spaced Data • One way to calculated derivatives of unequally spaced data is to determine a polynomial fit and take its derivative at a point. • As an example, using a second-order Lagrange polynomial to fit three points and taking its derivative yields:

Derivatives and Integrals for Data with Errors • A shortcoming of numerical differentiation is that it tends to amplify errors in data, whereas integration tends to smooth data errors. • One approach for taking derivatives of data with errors is to fit a smooth, differentiable function to the data and take the derivative of the function.

Numerical Differentiation with MATLAB • MATLAB has two built-in functions to help take derivatives, diff and gradient: • diff(x) – Returns the difference between adjacent elements in x • diff(y). /diff(x) – Returns the difference between adjacent values in y divided by the corresponding difference in adjacent values of x

Numerical Differentiation with MATLAB • fx = gradient(f, h) Determines the derivative of the data in f at each of the points. The program uses forward difference for the first point, backward difference for the last point, and centered difference for the interior points. h is the spacing between points; if omitted h=1. • The major advantage of gradient over diff is gradient’s result is the same size as the original data. • Gradient can also be used to find partial derivatives for matrices: [fx, fy] = gradient(f, h)

Visualization • MATLAB can generate contour plots of functions as well as vector fields. Assuming x and y represent a meshgrid of x and y values and z represents a function of x and y, – contour(x, y, z) can be used to generate a contour plot – [fx, fy]=gradient(z, h) can be used to generate partial derivatives and – quiver(x, y, fx, fy) can be used to generate vector fields

Chapter 20 Objectives • Understanding the meaning of local and global truncation errors and their relationship to step size for one-step methods for solving ODEs. • Knowing how to implement the following Runge-Kutta (RK) methods for a single ODE: – – Euler Heun Midpoint Fourth-Order RK • Knowing how to iterate the corrector of Heun’s method. • Knowing how to implement the following Runge-Kutta methods for systems of ODEs: – Euler – Fourth-order RK

Ordinary Differential Equations • Methods described here are for solving differential equations of the form: • The methods in this chapter are all one-step methods and have the general format: where is called an increment function, and is used to extrapolate from an old value yi to a new value yi+1.

Euler’s Method • The first derivative provides a direct estimate of the slope at ti: and the Euler method uses that estimate as the increment function:

Error Analysis for Euler’s Method • The numerical solution of ODEs involves two types of error: – Truncation errors, caused by the nature of the techniques employed – Roundoff errors, caused by the limited numbers of significant digits that can be retained • The total, or global truncation error can be further split into: – local truncation error that results from an application method in question over a single step, and – propagated truncation error that results from the approximations produced during previous steps.

Error Analysis for Euler’s Method • The local truncation error for Euler’s method is O(h 2) and proportional to the derivative of f(t, y) while the global truncation error is O(h). • This means: – The global error can be reduced by decreasing the step size, and – Euler’s method will provide error-free predictions if the underlying function is linear. • Euler’s method is conditionally stable, depending on the size of h.

MATLAB Code for Euler’s Method

Heun’s Method • One method to improve Euler’s method is to determine derivatives at the beginning and predicted ending of the interval and average them: • This process relies on making a prediction of the new value of y, then correcting it based on the slope calculated at that new value. • This predictor-corrector approach can be iterated to convergence:

Midpoint Method • Another improvement to Euler’s method is similar to Heun’s method, but predicts the slope at the midpoint of an interval rather than at the end: • This method has a local truncation error of O(h 3) and global error of O(h 2)

Runge-Kutta Methods • Runge-Kutta (RK) methods achieve the accuracy of a Taylor series approach without requiring the calculation of higher derivatives. • For RK methods, the increment function can be generally written as: where the a’s are constants and the k’s are where the p’s and q’s are constants.

Classical Fourth-Order Runge. Kutta Method • The most popular RK methods are fourthorder, and the most commonly used form is: where:

Systems of Equations • Many practical problems require the solution of a system of equations: • The solution of such a system requires that n initial conditions be known at the starting value of t.

Solution Methods • Single-equation methods can be used to solve systems of ODE’s as well; for example, Euler’s method can be used on systems of equations - the one-step method is applied for every equation at each step before proceeding to the next step. • Fourth-order Runge-Kutta methods can also be used, but care must be taken in calculating the k’s.

MATLAB RK 4 Code

Chapter 21 Objectives • Understanding how the Runge-Kutta Fehlberg methods use RK methods of different orders to provide error estimates that are used to adjust step size. • Familiarizing yourself with the built-in MATLAB function for solving ODEs. • Learning how to adjust options for MATLAB’s ODE solvers. • Learning how to pass parameters to MATLAB’s ODE solvers. • Understanding what is meant by stiffness and its implications for solving ODEs.

Adaptive Runge-Kutta Methods • The solutions to some ODE problems exhibit multiple time scales - for some parts of the solution the variable changes slowly, while for others there abrupt changes. • Constant step-size algorithms would have to apply a small stepsize to the entire computation, wasting many more calculations on regions of gradual change. • Adaptive algorithms, on the other hand, can change step-size depending on the region.

Approaches to Adaptive Methods • There are two primary approaches to incorporate adaptive step-size control: – Step halving - perform the one-step algorithm two different ways, once with a full step and once with two half-steps, and compare the results. – Embedded RK methods - perform two RK iterations of different orders and compare the results. This is the preferred method.

MATLAB Functions • MATLAB’s ode 23 function uses second- and thirdorder RK functions to solve the ODE and adjust step sizes. • MATLAB’s ode 45 function uses fourth- and fifthorder RK functions to solve the ODE and adjust step sizes. This is recommended as the first function to use to solve a problem. • MATLAB’s ode 113 function is a multistep solver useful for computationally intensive ODE functions.

Using ode Functions • The functions are generally called in the same way; ode 45 is used as an example: [t, y] = ode 45(odefun, tspan, y 0) – y: solution array, where each column represents one of the variables and each row corresponds to a time in the t vector – odefun: function returning a column vector of the righthand-sides of the ODEs – tspan: time over which to solve the system • If tspan has two entries, the results are reported for those times as well as several intermediate times based on the steps taken by the algorithm • If tspan has more than two entries, the results are reported only for those specific times – y 0: vector of initial values

Example - Predator-Prey • Solve: • • with y 1(0)=2 and y 2(0)=1 for 20 seconds predprey. m M-file: function yp = predprey(t, y) yp = [1. 2*y(1)-0. 6*y(1)*y(2); … -0. 8*y(2)+0. 3*y(1)*y(2)]; tspan = [0 20]; y 0 = [2, 1]; [t, y] = ode 45(@predprey, tspan, y 0); figure(1); plot(t, y); figure(2); plot(y(: , 1), y(: , 2));

ODE Solver Options • Options to ODE solvers may be passed as an optional fourth argument, and are generally created using the odeset function: options=odeset(‘par 1’, ‘val 1’, ‘par 2’, ‘val 2’, …) • Commonly used parameters are: – ‘Rel. Tol’: adjusts relative tolerance – ‘Abs. Tol’: adjusts absolute tolerance – ‘Initial. Step’: sets initial step size – ‘Max. Step’: sets maximum step size (default: one tenth of tspan interval)

Multistep Methods • Multistep methods are based on the insight that, once the computation has begun, valuable information from the previous points exists. • One example is the non-selfstarting Heun’s Method, which has the following predictor and corrector equations:

Stiffness • A stiff system is one involving rapidly changing components together with slowly changing ones. • An example of a single stiff ODE is: whose solution if y(0)=0 is:

MATLAB Functions for Stiff Systems • MATLAB has a number of built-in functions for solving stiff systems of ODEs, including ode 15 s, ode 23 t, and ode 23 tb. • The arguments for the stiff solvers are the same as those for previous solvers.