Remark: In many real-life mathematical models, one must utilize the integral or the derivative of a given function. Sometimes, these functions are easily differentiable or integrable. Many times they are not or doing so would be difficult. The question becomes determining how we can find good approximations to derivatives evaluated at given points and to definite integrals.
 

One idea: For any given function (or set of data points), determine a set of n+1 data points (x,y) and use the concepts from Chapter 3 to find a polynomial approximation to the original function. Usually, this polynomial approximation is very easy to differentiate and integrate. (Actually, this is not quite as true when using the Newton-type forms of the polynomial.) Use the resulting values from the polynomial as approximate values for the exact values.
 

Computing Assignment: Determine the Lagrange interpolating polynomial to f(x)=1/(1+x²) on the interval [-3,3] of degree n for various n. (Use the matrix approach by setting up the Vandermonde matrix, etc) Then, differentiate and integrate this polynomial exactly using formulas and compare with the exact values. Experiment with at least two other nontrivial functions for f(x) for which you can get an exact answer via calculus.
 

Remark: In the following, we will only consider determining derivatives at given points and definite integrals since each ot these is an actual number. If an interval approximation to the derivative function is desired or the indefinite integral is needed, the ideas which follow will need to be implemented over and over, once for each new value desired.
 

Derivatives over [a,b]:

1. Partition the interval [a,b] into x₀, x₁, ... ,x_n.
2. Use one of the following methods to determine an approximation m_k to the derivative for each x_k.
3. Connect these points (x_k,m_k ) with some interpolant from chapter 3.

1. Forward-Difference formula - page 140. (Difference quotient when h>0.) Has error O(h).
2. Backward-Difference formula - page 140. (Difference quotient when h<0.) Has error O(h).
3. Centered-Difference (Midpoint) formula - page 142. (Difference quotient centered around a.) Has error O(h²).
4. See various formulas on page 142, etc.
5. Second Derivative formula - page 146

Homework: page 146, #5, 6
 

Integrals over [a,b]:

1. Partition the interval [a,b] into x₀, x₁, ... ,x_n.
2. Use one of the following methods to determine an approximation I_k for the definite integral from x₀ to x_k .
• If I_k is known, then

I_k+1 = ,/`_{[a,x(k+1) ]} f(x) dx = F(x_k+1) - F(a)
      = F(x_k+1 ) - F(x_k) + F(x_k) - F(a)
      = F(x_k+1 ) - F(x_k) + I_k .
• Therefore, we can get the integrals data points recursively by integration over one more small interval at at time.
3. Connect these points (x_k,I_k ) with some interpolant from chapter 3.
4. Note that the answer might be off by an arbitrary constant.

Numerical Integration Methods:

From the definition of integration (Limits of Reimann sums), one perhaps would think that an approximation to the value of the definite integral could be found by replacing the infinite sum with a finite one. In general, numerical quadrature is approximating the definite integral with the finite sum:

S a_k f(x_k),

where the x_k are a reasonably fine finite partition of the interval [a,b].

1. Midpoint Rule: If the number of subintervals is even, on each pair of subintervals approximate the original function with a constant having the value of the function at the center of the subinterval (ie. a constant interpolant). Notice, this method does not require the function to be evaluated at either endpoint and so is a good choice to use if the integrand is singular at an endpoint.
2. Trapezoidal Rule: On each subinterval, approximate the "area under the curve" by using trapezoids. This is equivalent to approximating the actual function with a piecewise linear one connecting the function values at the endpoints of the little intervals. Lagrange gives the formula for linear interpolant.
3. Simpson's Rule: If the number of subintervals is even, on each pair of subintervals approximate the original function with a quadratic and find the area using this piecewise quadratic. Lagrange gives the formula for quadratic interpolant.
4. Gaussian Quadrature: Most accurate. Allow both coeffiecients a_k and x_k to be variable. Use table on page 110 for [a,b]=[-1,1].

Computing Assignments: Implement as separate functions the forward difference, backward difference and central difference methods for determining the derivative as well as the trapezoidal rule, Simpson's rule and the midpoint rule for determining the definite integral. Approximate the deriviative and integral of f(x) = 1/(1+x²) for several different values. Compare your answers with the known exact derivative f'(x) = -2x/(1+x ²)² and the (indefinite) integral tan^-1(x). Consider other functions also and compare your approximation with the exact value via calculus.

Ex: Solving y'=4y gives y(t)=Ce^4t , where C is any constant. To eliminate the arbitrary constant C, we would need an initial condition for y(a) given. For example, y(0)=3 implies y(t)=3e ^4t.
 

General Problem: Solve for y(t) on the interval a < t < b,

(IVP) y'=f(t,y(t))

y(a) = y₀ , given.
 

Remark: Elementary Ordinary Differential Equation (ODE) theory gives mathematical discussions on the existence and uniqueness of exact solutions to (IVP). (See Well-Posed Condition on page 151.) We will always assume that we are given problems for which solutions exist and are unique.
 

Remark: If (IVP) cannot be solved by analytical methods, we will need to find a discrete solution by finding an approximation to the solution at a finite number of points. Once we have these points, we can use approximation theory from Chapter 3 to find an approximate continuous solution.
 

Notation: Given interval endpoints a and b; number of data points desired n+1; and the function f(t,y) from the D.E....

Let h=(b-a)/n, be the mesh size with equally spaced mesh points given by t_j=a+j h; w_j= algorithm generated approximation for y(t_j); and set f_j=f(t _j,w_j), for j=0,1,2,...,n.
 

Euler's Method:

w₀ = y ₀,

w_j+1 = fw_j+h*f_j , for j=0,1,2,...n.

Derivations of Euler's Method:

(1) Take y'(t_j) = f(t_j,y(t _j)) and replace y(t_j) by a forward difference formula.

(2) Use Taylor's theorem to expand y(t_j+1 ) about t_j and throw away nonlinear terms. Easily gives local error. Global Error, pg 154.

(3) Integrate from t_j to t_j+1 and approximate f(u,y(u)) under the integral by a constant evaluated at u=t _j.

(4) Graphically, find the equation of the tangent line at t=t_j and evaluate it at t=t_j+1 to approximate y(t_j+1).

But m=y'=f(t_j,y(t_j)) gives y-w_j = f_j(t - t_j). Plugging in t=t_j+1 gives the method.
 

Assignment: Program Euler's Method in for solving the Differential Equations above. Represent the numerical solution graphically and compare it graphically with the known exact solutions. Give output corresponding to several different values for the step size h.

Remark: We would like to compare methods for solving (IVP) to see which are "better" to use. To do this, consider the difference equation written in the form w_k+1 = w_j + h*f(t_k,w_k), where f is "all the rest" of the method. Define the local truncation error of a method to be t(h) = (y(t _k+1) - y(t_k))/h + f(t_k,y_k). The order of the method is the lowest power of h remaining in t(h) after expanding all terms about t=t_k and simplifying.
 

Result: Euler's method is order 1.
 

Taylor's Methods: page 156

To develop methods of order higher than one, use the Taylor's exansion for y(t_k+h) about t_k and note that the derivatives of y can be obtained to arbitrarly high degree by using the D.E. y'=f(t,y), y''=f'(t,y), y'''=f''(t,y), etc, where f(t,y(t)) is a given function.
 

Homework: page 159, #3, 4, 7.
 

Runge Kutta Methods:

Instead of evaluating f(t,y) only at "nodal" points, allow this function to be evaluated at points between nodes. The extra degrees of freedom give higher order methods which don't involve the derivative calculations needed in Taylor's methods.

(a) Midpoint Method - page 163

(b) Modified Euler Method - page 163

(c) Heun's Method - page 163

(d) Classical Runge-Kutta Method - page 164 (This is the method most often used in practice.)
 

Homework: page 166, #4, 11 (Feel free to use Software.)
 

Multi-Step Methods: page 168

Instead of only using information at time t to get the approximation at time t+h, use the computational values from previous time steps also. This enables one to get methods of high accuracy without the complexity of using the Runge-Kutta methods. The tradeoff is that this type of method needs more that one starting value and these must be generated by a single-step method.

-Adam's-Bashforth Methods - page 169, explicitly gives w_k+1.

-Adam's-Moulton Methods - page 170, implicitly gives w_k+1.

Often, these multi-step methods are used together with an Adam's-Bashforth method used to give a first approximation to w _k+1 and then an Adam's-Moulton method used to improve this value.
 

Computing Assignment:

Implement both the Midpoint and Classical Runge-Kutta methods.
 

Optional Programming Assignment:

Implement the Adam's Fourth-Order Predictor-Corrector Algorithm on page 266.

Discuss:
 

1. System of equations
2. Matrix notation and the operations on matrices.
3. Solving Diagonal Systems and Triangular Systems.
• Forward Solving - O(n²)
• Backward Solving - O(n²)
4. Gaussian Elimination.
5. Factorization of Matrices and using this to obtain the determinant in O(n³) time (vs. n! time with cofactor expansion)
6. Other factorizations: CROUT, DOOLITTLE, CHOLESKY
7. Triangular Solvers.

(a) Determination of A^-1 using the factorization. Notice, if z_k is the kth column of A^-1 and e_k is the column vector with all zeros except for a one in the kth position, we can write the matrix equation A*A^-1 = I as n matrix equations A*z_k = e_k. Factoring the matrix A into A=LU once, yields LUz_k=e_k, and notice each of these can be solved in O(n²) time.

(b) Consider a tensor product interpolating surface x(u,v) = S S c_i,j A_i(u) B_j(v), where A _i(u) and B_j(v) are basis functions. Suppose we are given data points to interpolate x_i,j. Plugging these points into the above tensor product yields a system of equations X=ACB, where X, A, C, B are matrices. The formal solution of this is C=A^-1XB^-1 , which requires O(n⁶) operations. Using a factorization process carefully yields only O(n⁴) operations - a significant savings for even small systems.
 
 

Check for correctness, legibility, test cases, presentation in MathCad, look for programming extra's or extra effort.
 

_____ Bernstein
 
 
 
 
 

_____ Bisection
 
 
 
 
 

_____ Newtons
 
 
 
 
 

_____ Secant
 
 
 
 
 

_____ Newton's Divided Differences
 
 
 
 
 

_____ MathCad Splines
 
 
 
 
 

_____ Lagrange Polynomial Derivatives and Integrals
 
 
 
 
 

_____ Difference Formulas
 
 
 
 
 

_____ Quadrature Formulas
 
 
 
 
 

_____ Euler's Method
 
 
 
 
 

_____ 2nd Order Taylor's Method
 
 
 
 
 

_____ Midpoint Runge Kutta
 
 
 
 
 

_____ Classical Runge Kutta