Introductory Methods of Numerical Analysis

Runge-Kutta-Fehlberg Method

2010-02-02T04:13:00.000-08:00

One way to guarantee accuracy in the solution of an I.V.P. is to solve the problem twice using step sizes h and and compare answers at the mesh points corresponding to the larger step size. But this requires a significant amount of computation for the smaller step size and must be repeated if it is determined that the agreement is not good enough. The Runge-Kutta-Fehlberg method (denoted RKF45) is one way to try to resolve this problem. It has a procedure to determine if the proper step size h is being used. At each step, two different approximations for the solution are made and compared. If the two answers are in close agreement, the approximation is accepted. If the two answers do not agree to a specified accuracy, the step size is reduced. If the answers agree to more significant digits than required, the step size is increased.

Each Runge-Kutta-Fehlberg step requires the use of the following six values:

Then an approximation to the solution of the I.V.P. is made using a Runge-Kutta method of order 4:

And a better value for the solution is determined using a Runge-Kutta method of order 5:

The optimal step size sh can be determined by multiplying the scalar s times the current step size h. The scalar s is

where is the specified error control tolerance.

Adams-Bashforth-Moulton Method

2010-02-02T04:07:00.000-08:00

The methods of Euler, Heun, Taylor and Runge-Kutta are called single-step methods because they use only the information from one previous point to compute the successive point, that is, only the initial point is used to compute and in general is needed to compute . After several points have been found it is feasible to use several prior points in the calculation. The Adams-Bashforth-Moulton method uses in the calculation of . This method is not self-starting; four initial points , , , and must be given in advance in order to generate the points .

A desirable feature of a multistep method is that the local truncation error (L. T. E.) can be determined and a correction term can be included, which improves the accuracy of the answer at each step. Also, it is possible to determine if the step size is small enough to obtain an accurate value for , yet large enough so that unnecessary and time-consuming calculations are eliminated. If the code for the subroutine is fine-tuned, then the combination of a predictor and corrector requires only two function evaluations of f(t,y) per step.

Adams-Bashforth-MoultonMethod:

Assume that f(t,y) is continuous and satisfies a Lipschits condition in the variable y, and consider the I. V. P. (initial value problem)

with , over the interval .

The Adams-Bashforth-Moulton method uses the formulas , and

the predictor , and

the corrector for

as an approximate solution to the differential equation using the discrete set of points .

Remark: The Adams-Bashforth-Moulton method is not a self-starting method. Three additional starting values must be given. They are usually computed using the Runge-Kutta method.

Precision of Adams-Bashforth-MoultonMethod: Assume that is the solution to the I.V.P. with . If and is the sequence of approximations generated by Adams-Bashforth-Moulton method, then at each step, the local truncation error is of the order , and the overall global truncation error is of the order

, for .

The error at the right end of the interval is called the final global error

.

Adams-Bashforth-Moulton Method:

To approximate the solution of the initial value problem with over at a discrete set of points using the formulas:

use the predictor

and the corrector for .

Milne-Simpson's Method

2010-02-02T03:55:00.000-08:00

The methods of Euler, Heun, Taylor and Runge-Kutta are called single-step methods because they use only the information from one previous point to compute the successive point, that is, only the initial point is used to compute and in general is needed to compute . After several points have been found it is feasible to use several prior points in the calculation. The Milne-Simpson method uses in the calculation of . This method is not self-starting; four initial points , , , and must be given in advance in order to generate the points .

A desirable feature of a multistep method is that the local truncation error (L. T. E.) can be determined and a correction term can be included, which improves the accuracy of the answer at each step. Also, it is possible to determine if the step size is small enough to obtain an accurate value for , yet large enough so that unnecessary and time-consuming calculations are eliminated. If the code for the subroutine is fine-tuned, then the combination of a predictor and corrector requires only two function evaluations of f(t,y) per step.

Theorem (Milne-Simpson's Method:

Assume that f(t,y) is continuous and satisfies a Lipschits condition in the variable y, and consider the I. V. P. (initial value problem)

with , over the interval .

The Milne-Simpson method uses the formulas , and

the predictor , and

the corrector for

as an approximate solution to the differential equation using the discrete set of points .

Remark: The Milne-Simpson method is not a self-starting method. Three additional starting values must be given. They are usually computed using the Runge-Kutta method.

Theorem (Precision of the Milne-Simpson Method:

Assume that is the solution to the I.V.P. with . If and is the sequence of approximations generated by Milne-Simpson method, then at each step, the local truncation error is of the order , and the overall global truncation error is of the order

, for .

The error at the right end of the interval is called the final global error

.

Shooting Methods for ODE's

2010-02-02T03:50:00.000-08:00

Finding the solution of a linear boundary problem is assisted by the linear structure of the equation and the use of two special initial value problems. Suppose that u(t) is the unique solution to the I.V.P.

with .

Furthermore, suppose that v(t) is the unique solution to the I.V.P.

with .

Then the linear combination

.

is a solution to with .

Program (Linear Shooting Method):

To approximate the solution of the boundary value problem

with

over the interval [a,b] by using the Runge-Kutta method of order n=4.

The method involves solving a two systems of equations over . First solve

with ,
and .

Then solve

with ,
and .

Finally, the desired solution x(t) is the linear combination

.

The subroutine Runge2D will be used to construct the two solutions , and .

Finite Difference Method for ODE's

2010-02-02T03:43:00.000-08:00

Methods involving difference quotient approximations for derivatives can be used for solving certain second-order boundary value problems. Consider the linear equation

(1)

over [a,b] with . Form a partition of [a, b] using the points , where and for . The central-difference formulas discussed in Chapter 6 are used to approximate the derivatives

(2)

and

(3)

Use the notation for the terms on the right side of (2) and (3) and drop the two terms . Also, use the notations , , and this produces the difference equation

which is used to compute numerical approximations to the differential equation (1). This is carried out by multiplying each side by and then collecting terms involving and arranging them in a system of linear equations:

for , where and . This system has the familiar tridiagonal form.

Galerkin's Method

2010-02-02T03:32:00.000-08:00

One of the most important weighted residual methods was invented by the Russian mathematician Boris Grigoryevich Galerkin (February 20, 1871 - July 12, 1945). Galerkin's method selects the weight function functions in a special way: they are chosen from the basis functions, i.e. . It is required that the following equations hold true

(6) for .

To apply the method, all we need to do is solve these equations for the coefficients .

Galerkin's Method for solving an I. V. P.

Suppose we wish to solve the initial value problem

(i) ,
with
over the interval .

We use the trial function

(ii) .

There are equations to solve for , i.e.

(iii) for .

Remark: For the solution of an I. V. P. we choose .

Galerkin's Method for solving an a B. V. P.

Suppose we wish to solve a boundary value problem over the interval ,

(I) ,
with

We define and use the trial function

(II) .

There are equations to solve for , i.e.

(III) for .

Remark: The functions must all be chosen with the boundary properties

and for .

Lotka-Volterra Model

2010-02-02T03:24:00.000-08:00

Lotka-Volterra Equations:
The "Lotka-Volterra equations" refer to two coupled differential equations

There is one critical point which occurs when and it is .

The Runge-Kutta method is used to numerically solve O.D.E.'s over an interval .

Pendulum

2010-02-02T03:13:00.000-08:00

Nonlinear Pendulum:
A simple pendulum consists of a point mass m attached to a rod of negligible weight. The torque is

(1) ,

where denotes the the angle of the rod measured downward from a vertical axis. The moment of inertia for the point mass is where l is the length of the rod. The torque can also be expressed as , where is the angular acceleration, using Newton's second law, and the second derivative, this can be written as

(2) .

Equating (1) and (2) results in the nonlinear D. E.

(3) .

Linear Pendulum:
Introductory courses discuss the pendulum with small oscillations as an example of a simple harmonic oscillator. If the angle of oscillation is small, use the approximation in equation (3) and obtain the familiar linear D. E. for simple harmonic motion:

(4) ,

Using the substitution , the solution to (4) is known to be

(5) ,

which has period . When the solution (5) is written with a phase shift, it becomes

(6) .

Projectile Motion

2010-02-02T03:01:00.000-08:00

In calculus, a model for projectile motion with no friction is considered, and a "parabolic trajectory" is obtained. If the initial velocity is and is the initial angle to the horizontal, then the parametric equations for the horizontal and vertical components of the position vector are

(1) ,
and
(2) .

Solve equation (1) for t and get , then replace this value of t in equation (2) and the result is

,

which is an equation of a parabola.

The time required to reach the maximum height is found by solving :

,
yields
,

and the maximum height is

.

The time till impact is found by solving , which yields , and for this model, . The range is found by calculating :

.

For a fixed initial velocity , the range is a function of , and is maximum when .

Numerical solution of second order D. E.'s

This module illustrates numerical solutions of a second order differential equation. First, we consider the special case where the projectile is fired vertically along the y-axis and has no horizontal motion, i. e. x(t)= 0. The effect of changing the amount of air drag or air resistance is investigated. It is known that the drag force acting on an object which moves very slowly through a viscous fluid is directly proportional to the velocity of that object. However, there are examples, such as Millikan's oil drop experiment, when the drag force is proportional to the square of the velocity. Further investigations into the situation could involve the Reynolds number.

Math-Models (Projectile Motion I):

The following mathematical models are are considered.
(i). No air resistance , and
.

(ii). Air resistance proportional to velocity , and
.

(iii). Air resistance proportional to the square of the velocity , for the ascent, and , for the descent, and

(iv). Air resistance proportional to the power of the velocity

Lorenz Attractor

2010-02-02T02:30:00.000-08:00

The Lorenz attractor is an attractor that arises in a simplified system of equations describing the two-dimensional flow of fluid of uniform depth , with an imposed temperature difference , under gravity , with buoyancy , thermal diffusivity , and kinematic viscosity . The full equations are

			(1)
			(2)

Here, is a stream function, defined such that the velocity components of the fluid motion are

			(3)
			(4)

(Tabor 1989, p. 205).

In the early 1960s, Lorenz accidentally discovered the chaotic behavior of this system when he found that, for a simplified system, periodic solutions of the form

(5)

(6)

grew for Rayleigh numbers larger than the critical value, Ra_c" border="0" height="16" width="58">. Furthermore, vastly different results were obtained for very small changes in the initial values, representing one of the earliest discoveries of the so-called butterfly effect.

Lorenz included the terms

			(7)
			(8)
			(9)

where is proportional to convective intensity, to the temperature difference between descending and ascending currents, and to the difference in vertical temperature profile from linearity in his system of equations. From these, he obtained the simplified equations

			(10)
			(11)
			(12)

now known as the Lorenz equations. Here, , , , and

			(13)
			(14)
			(15)

where is the Prandtl number, Ra is the Rayleigh number, is the critical Rayleigh number, and is a geometric factor (Tabor 1989, p. 206). Lorenz took and .

The Lorenz attractor has a correlation exponent of and capacity dimension (Grassberger and Procaccia 1983). For more details, see Lichtenberg and Lieberman (1983, p. 65) and Tabor (1989, p. 204). As one of his list of challenging problems for mathematics (Smale's problems), Smale (1998, 2000) posed the open question of whether the Lorenz attractor is a strange attractor. This question was answered in the affirmative by Tucker (2002), whose technical proof makes use of a combination of normal form theory and validated interval arithmetic.

The critical points at (0, 0, 0) correspond to no convection, and the critical points at

(16)

and

(17)

correspond to steady convection. This pair is stable only if

(18)

which can hold only for positive if b+1" border="0" height="14" width="54">.

The image above shows a Lorenz attractor laser-etched into glass by digital sculptor Bathsheba Grossman

The Lorenz attractor is a set of differential equations which are popular in the field of Chaos. The equations describe the flow of fluid in a box which is heated along the bottom. This model was intended to simulate medium-scale atmospheric convection. Lorenz simplified some of the Navier-Stokes equations in the area of fluid dynamics and obtained three ordinary differential equations

,

.

The parameter p is the Prandtl number, is the quotient of the Rayleigh number and critical Rayleigh number and b is a geometric factor. Lorenz is attributed to using the values .

There are three critical points (0,0,0) corresponds to no convection, and the two points

and correspond to steady convection.

The latter two points are to be stable, only if the following equation holds

.

Van Der Pol System

2010-02-02T02:18:00.000-08:00

The van der Pol equation is an ordinary differential equation that can be derived from the Rayleigh differential equation by differentiating and setting . It is an equation describing self-sustaining oscillations in which energy is fed into small oscillations and removed from large oscillations. This equation arises in the study of circuits containing vacuum tubes and is given by

If , the equation reduces to the equation of simple harmonic motion

The van der Pol equation is

,
where
is a constant.

When the equation reduces to , and has the familiar solution . Usually the term in equation (1) should be regarded as friction or resistance, and this is the case when the coefficient is positive. However, if the coefficient is negative then we have the case of "negative resistance." In the age of "vacuum tube" radios, the "tetrode vacuum tube" (cathode, grid, plate), was used for a power amplifier and was known to exhibit "negative resistance." The mathematics is amazing too, and van der Pol, Balthasar (1889-1959) is credited with developing equation (1). The solution curves exhibits orbital stability. The van der Pol equation can be written as a second order system

,
and
.

Any convenient numerical differential equation solver such as the Runge-Kutta method can be used to compute the solutions.

Harvesting Model

2010-02-02T02:08:00.000-08:00

Harvesting a Logistic Population:

When the harvesting term -k is incorporated into into bounded population model we have

.

There are three solution forms for this differential equation, and they correspond to the nature of the stationary solutions ( x(t) = c).

Definition(Stationary Points):

The stationary points of the D. E. are solutions where and are the roots of the characteristic equation

.

The roots are known to be , and the stationary solutions are .

Remark:

Since x(t) is a real function, there are no stationary solutions when .

Case (i) If there is one stationary solution :

When , the differential equation has the form and the solution is

.

The solution with the initial condition is

.

If then .

If then function x(t) has a vertical asymptote at

and the population x(t) becomes extinct at some time (where ), i. e.

.

Case (ii) If there are two stationary solutions and :

When , the differential equation has the form and the solution is

.

The two real roots of the characteristic equation , are .

The solution with the initial condition

If then .

If then the population x(t) becomes extinct at some time , i. e. .

Case (iii) If there are no stationary solutions:

When , the differential equation has the form and the solution is

.

The solution with the initial condition is

The function x(t) has a vertical asymptote at so the population x(t) becomes extinct at some time (where .), i.e.

.

Frobenius Series Solution

2010-02-02T01:58:00.000-08:00

Consider the second order linear differential equation

(1) .

Rewrite this equation in the form , then use the substitutions and and rewrite the differential equation (1) in the form

(2) .

Definition (Analytic):

The functions and are analytic at if they have Taylor series expansions with radius of convergence and , respectively. That is

which converges for
and
which converges for

Definition (Ordinary Point):

If the functions and are analytic at , then the point is called an ordinary point of the differential equation

.

Otherwise, the point is called a singular point.

Definition (Regular Singular Point):

Assume that is a singular point of (1) and that and are analytic at .

They will have Maclaurin series expansions with radius of convergence and , respectively. That is

which converges for
and
which converges for

Then the point is called a regular singular point of the differential equation (1).

Method of Frobenius:

This method is attributed to the german mathemematican Ferdinand Georg Frobenius (1849-1917 ). Assume that is regular singular point of the differential equation

.

A Frobenius series (generalized Laurent series) of the form

can be used to solve the differential equation. The parameter must be chosen so that when the series is substituted into the D.E. the coefficient of the smallest power of is zero. This is called the indicial equation. Next, a recursive equation for the coefficients is obtained by setting the coefficient of equal to zero. Caveat: There are some instances when only one Frobenius solution can be constructed.

Definition (Indicial Equation):

The parameter in the Frobenius series is a root of the indicial equation

.

Assuming that the singular point is , we can calculate as follows:

and

The Recursive Formulas:

For each root of the indicial equation, recursive formulas are used to calculate the unknown coefficients . This is custom work because a numerical value for is easier use.

Picard Iteration

2010-02-02T01:47:00.000-08:00

The term "Picard iteration" occurs two places in undergraduate mathematics. In numerical analysis it is used when discussing fixed point iteration for finding a numerical approximation to the equation . In differential equations, Picard iteration is a constructive procedure for establishing the existence of a solution to a differential equation that passes through the point .

The first type of Picard iteration uses computations to generate a "sequence of numbers" which converges to a solution. We will not present this application, but mention that it involves the traditional role of computer software as a "number cruncher."

The goal of this article is to illustrate the second application of Picard iteration; i. e. how to use a computer to efficiently generate a "sequence of functions" which converges to a solution. We will see that computer software can perform the more sophisticated task of "symbol cruncher

Most differential equations texts give a proof for the existence and uniqueness of the solution to a first order differential equation. Then exercises are given for performing the laborious details involved in the method of successive approximations. The concept seems straightforward, just repeated integration, but students get bogged down with the details. Now computers can do all the drudgery and we can get a better grasp on how the process works.

Theorem 1 (Existence Theorem):
If both are continuous on the rectangle and , then there exists a unique solution to the initial value problem (I.V.P.)

(1)

for all values of x in some (smaller) interval contained in .

Picard's Method for D.E.'s

The method of successive approximations uses the equivalent integral equation for (1) and an iterative method for constructing approximations to the solution. This is a traditional way to prove (1) and appears in most all differential equations textbooks. It is attributed to the French mathematician Charles Emile Picard (1856-1941).

Theorem 2 (Successive Approximations - Picard Iteration):

The solution to the I.V.P in (1) is found by constructing recursively a sequence of functions

, and
(2)
.

Then the solution to (1) is given by the limit:

(3) .

Spring-Mass Systems

2010-02-01T10:07:00.000-08:00

Consider the system of two masses and two springs with no external force. Visualize a wall on the left and to the right a spring , a mass, a spring and another mass. Assume that the spring constants are . See Figure 1 below.

Figure 1. Coupled masses with spring attached to the wall at the left.

Assume that the masses slide on a frictionless surface and that the functions denote the displacement from static equilibrium of the masses , respectively. It can be shown by using Newton's second law and Hooke's law that the system of D. E.'s for is

Remark:

The eigenfrequencies can be obtained by taking the square root of the eigenvalues of the matrix .

Consider the system of two masses and three springs with no external force. Visualize a wall on the left and to the right a spring , a mass, a spring, a mass, a spring and another wall. Assume that the spring constants are . See Figure 2 below.

Figure 2. Coupled masses with springs attached to walls at the left and right.

Assume that the masses slide on a frictionless surface and that the functions denote the displacement from static equilibrium of the masses , respectively. It can be shown by using Newton's second law and Hooke's law that the system of D. E.'s for is

Remark:

The eigenfrequencies can be obtained by taking the square root of the eigenvalues of the matrix .

Eigen Frequencies:

Consider three springs in parallel, with two of the springs having spring constant k and attached to two walls on either end, and the third spring of spring constant k placed between two equal masses m. To solve for the motion of the masses using the normal formalism, equate forces


			(1)

			(2)

Writing (1) and (2) in matrix form gives

(3)

To find the equations of motion from energy considerations, note that the kinetic energy is defined by

(4)

so the rate of change of K is

(5)

Similarly, the potential energy is defined by


			(6)

so the rate of change of U is



			(7)

But energy is conserved, so

(8)

and

(9)

Now, this equation must hold for arbitrary and , so each piece must vanish separately ("separation of variables "), yielding the coupled equations (3). The eigenmodes of the system follow from (3). Looking for a harmonic solution using the trial solution ,


			(10)

(11)

This has solutions when the determinant is 0, so

(12)

(13)

(14)

To find the eigenvectors, plug back in. For ,

(15)

so and the first eigenvalue and its associated eigenvector are

(16)

This corresponds to masses moving in opposite directions. For ,

(17)

(18)

This corresponds to masses moving in the same direction.

If the middle spring has the same spring constant as those on either side, then , and the eigenfrequencies are

			(19)
			(20)

Crank-Nicolson Method

2010-02-01T09:55:00.000-08:00

An implicit scheme, invented b John Crank (1916-) and Phyllis Nicolson (1917-1968), is based on numerical approximations for solutions at the point that lies between the rows in the grid. Specifically, the approximation used for is obtained from the central-difference formula,

.

Elliptic Partial Differential Equations

2010-02-01T09:35:00.000-08:00

As examples of elliptic partial differential equations, we consider the Laplace equation, Poisson equation, and Helmholtz equation. Recall that the Laplacian of the function u(x,y) is

.

With this notation, we can write the Laplace, Poisson, and Helmholtz equations in the following forms:

It is often the case that the boundary values for the function u(x,y) are known at all points on the sides of a rectangular region R in the plane. In this case, each of these equations can be solved by the numerical technique known as the finite-difference method.

Eigenvalues and Eigenvectors

2010-02-01T09:24:00.000-08:00

We will now review some ideas from linear algebra. Proofs of the theorems are either left as exercises or can be found in any standard text on linear algebra. We know how to solve n linear equations in n unknowns. It was assumed that the determinant of the matrix was nonzero and hence that the solution was unique. In the case of a homogeneous system AX = 0, if , the unique solution is the trivial solution X = 0. If , there exist nontrivial solutions to AX = 0. Suppose that , and consider solutions to the homogeneous linear
system

A homogeneous system of equations always has the trivial solution . Gaussian elimination can be used to obtain the reduced row echelon form which will be used to form a set of relationships between the variables, and a non-trivial solution.

Definition (Linearly Independent):

The vectors are said to be linearly independent if the equation

implies that . If the vectors are not linearly independent they are said to be linearly dependent.

Two vectors in are linearly independent if and only if they are not parallel. Three vectors in are linearly independent if and only if they do not lie in the same plane.

Definition (Linearly Dependent):

The vectors are said to be linearly dependent if there exists a set of numbers not all zero, such that

.

Theorem:

The vectors are linearly dependent if and only if at least one of them is a linear combination of the others.

A desirable feature for a vector space is the ability to express each vector as s linear combination of vectors chosen from a small subset of vectors. This motivates the next definition.

Definition (Basis):

Suppose that is a set of m vectors in . The set S i s called a basis for if for every vector there exists a unique set
of scalars so that X can be expressed as the linear combination

Theorem:

In , any set of n linearly independent vectors forms a basis of . Each vector is uniquely expressed as a linear combination of the basis vectors.

Theorem:

Let be vectors in .

(i) If m>n, then the vectors are linearly independent.

(ii) If m=n, then the vectors are linearly dependent if and only if , where .

Applications of mathematics sometimes encounter the following questions: What are the singularities of , where is a parameter? What is the behavior of the sequence of vectors ? What are the geometric features of a linear transformation? Solutions for problems in many different disciplines, such as economics, engineering, and physics, can involve ideas related to these equations. The theory of eigenvalues and eigenvectors is powerful enough to help solve these otherwise intractable problems.

Let A be a square matrix of dimension n × n and let X be a vector of dimension n. The product Y = AX can be viewed as a linear transformation from n-dimensional space into itself. We want to find scalars for which there exists a nonzero vector X such that

(1) ;

that is, the linear transformation T(X) = AX maps X onto the multiple . When this occurs, we call X an eigenvector that corresponds to the eigenvalue , and together they form the eigenpair for A. In general, the scalar and vector X can involve complex numbers. For simplicity, most of our illustrations will involve real calculations. However, the techniques are easily extended to the complex case. The n × n identity matrix I can be used to write equation (1) in the form

(2) .

The significance of equation (2) is that the product of the matrix and the nonzero vector X is the zero vector! The theorem of homogeneous linear system says that (2) has nontrivial solutions if and only if the matrix is singular, that is,

(3) .

This determinant can be written in the form

(4)

Definition (Characteristic Polynomial):

When the determinant in (4) is expanded, it becomes a polynomial of degree n, which is called the characteristic polynomial

(5)

Exploration For p()

There exist exactly n roots (not necessarily distinct) of a polynomial of degree n. Each root can be substituted into equation (3) to obtain an underdetermined system of equations that has a corresponding nontrivial solution vector X. If is real, a real eigenvector X can be constructed. For emphasis, we state the following definitions.

Definition (Eigenvalue):

If A is and n × n real matrix, then its n eigenvalues are the real and complex roots of the characteristic polynomial

.

Definition (Eigenvector):

If is an eigenvalue of A and the nonzero vector V has the property that

then V is called an eigenvector of A corresponding to the eigenvalue . Together, this eigenvalue and eigenvector V is called an eigenpair .

The characteristic polynomial can be factored in the form

where is called the multiplicity of the eigenvalue . The sum of the multiplicities of all eigenvalues is n; that is,

.

The next three results concern the existence of eigenvectors.

Theorem (Corresponding Eigenvectors):

Suppose that A is and n × n square matrix.
(a) For each distinct eigenvalue there exists at least one eigenvector V corresponding to .
(b) If has multiplicity r, then there exist at most r linearly independent eigenvectors that correspond to .

Theorem (Linearly Independent Eigenvectors):

Suppose that A is and n × n square matrix. If the eigenvalues are distinct and are the k eigenpairs, then is a set of k linearly independent vectors.

Theorem (Complete Set of Eigenvectors):

Suppose that A is and n × n square matrix. If the eigenvalues of A are all distinct, then there exist n nearly independent eigenvectors .

Finding eigenpairs by hand computations is usually done in the following manner. The eigenvalue of multiplicity r is substituted into the equation

.

Then Gaussian elimination can be performed to obtain the row reduced echelon form, which will involve n-k equations in n unknowns, where . Hence there are k free variables to choose. The free variables can be selected in a judicious manner to produce k linearly independent solution vectors that correspond to .

Power Method

2010-02-01T08:53:00.000-08:00

We now describe the power method for computing the dominant eigenpair. Its extension to the inverse power method is practical for finding any eigenvalue provided that a good initial approximation is known. Some schemes for finding eigenvalues use other methods that converge fast, but have limited precision. The inverse power method is then invoked to refine the numerical values and gain full precision. To discuss the situation, we will need the following definitions.

If is an eigenvalue of A that is larger in absolute value than any other eigenvalue, it is called the dominant eigenvalue. An eigenvector corresponding to is called a dominant eigenvector.

An eigenvector V is said to be normalized if the coordinate of largest magnitude is equal to unity (i.e., the largest coordinate in the vector V is the number 1).

Remark:

It is easy to normalize an eigenvector by forming a new vector where and .

Theorem (Power Method):

Assume that the n×n matrix A has n distinct eigenvalues and that they are ordered in decreasing magnitude; that is, . If is chosen appropriately, then the sequences and generated recursively by

and

where and , will converge to the dominant eigenvector and eigenvalue , respectively. That is,

and .

Remark:

If is an eigenvector and , then some other starting vector must be chosen.

Speed of Convergence

In the iteration in the theorem uses the equation

,

and the coefficient of that is used to form goes to zero in proportion to . Hence, the speed of convergence of to is governed by the terms . Consequently, the rate of convergence is linear. Similarly, the convergence of the sequence of constants to is linear. The Aitken method can be used for any linearly convergent sequence to form a new sequence,

,

that converges faster. The Aitken can be adapted to speed up the convergence of the power method.

Shifted-Inverse Power Method

We will now discuss the shifted inverse power method. It requires a good starting approximation for an eigenvalue, and then iteration is used to obtain a precise solution. Other procedures such as the QM and Givens’ method are used first to obtain the starting approximations. Cases involving complex eigenvalues, multiple eigenvalues, or the presence of two eigenvalues with the same magnitude or approximately the same
magnitude will cause computational difficulties and require more advanced methods. Our illustrations will focus on the case where the eigenvalues are distinct. The shifted inverse power method is based on the following three results (the proofs are left as exercises).

Theorem (Shifting Eigenvalues):

Suppose that ,V is an eigenpair of A. If is any constant, then ,V is an eigenpair of the matrix .

Theorem (Inverse Eigenvalues):

Suppose that ,V is an eigenpair of A. If , then ,V is an eigenpair of the matrix .

Theorem (Shifted-Inverse Eigenvalues):

Suppose that ,V is an eigenpair of A. If , then ,V is an eigenpair of the matrix .

Theorem (Shifted-Inverse Power Method):

Assume that the n×n matrix A has distinct eigenvalues and consider the eigenvalue . Then a constant can be chosen so that is the dominant eigenvalue of . Furthermore, if is chosen appropriately, then the sequences and generated recursively by

and

where and , will converge to the dominant eigenpair , of the matrix . Finally, the corresponding eigenvalue for the matrix A is given by the calculation

Remark. For practical implementations of this Theorem, a linear system solver is used to compute in each step by solving the linear system .

Jacobi method

2010-02-01T08:46:00.000-08:00

Jacobi’s method is an easily understood algorithm for finding all eigenpairs for a symmetric matrix. It is a reliable method that produces uniformly accurate answers for the results. For matrices of order up to 10×10, the algorithm is competitive with more sophisticated ones. If speed is not a major consideration, it is quite acceptable for matrices up to order 20×20. A solution is guaranteed for all real symmetric matrices when Jacobi’s method is used. This limitation is not severe since many practical problems of applied mathematics and engineering involve symmetric matrices. From a theoretical viewpoint, the method embodies techniques that are found in more sophisticated algorithms. For instructive purposes, it is worthwhile to investigate the details of Jacobi’s method.

Jacobi Series of Transformations

Start with the real symmetric matrix . Then construct the sequence of orthogonal matrices as follows:

and
for j = 1, 2, ... .

It is possible to construct the sequence so that

.

In practice we will stop when the off-diagonal elements are close to zero. Then we will have

.

Current research by James W. Demmel and Kresimir Veselic (1992) indicate that Jacobi's method is more accurate than QR. You can check out their research by following the link in the list of internet resources. The abstract for their research follows below.

We show that Jacobi's method (with a proper stopping criterion) computes small eigenvalues of symmetric positive definite matrices with a uniformly better relative accuracy bound than QR, divide and conquer, traditional bisection, or any algorithm which first involves tridiagonalizing the matrix. In fact, modulo an assumption based on extensive numerical tests, we show that Jacobi's method is optimally accurate in the following sense: if the matrix is such that small relative errors in its entries cause small relative errors in its eigenvalues, Jacobi will compute them with nearly this accuracy. In other words, as long as the initial matrix has small relative errors in each component, even using infinite precision will not improve on Jacobi (modulo factors of dimensionality). ...

Householder's Method

2010-02-01T08:30:00.000-08:00

If and are vectors with the same norm, there exists an orthogonal symmetric matrix such that

,
where

and
.

Since is both orthogonal and symmetric, it follows that

.

Corollary ( Householder Matrix)

Let be an matrix, and any vector. If is an integer with , we can construct a vector and matrix so that

(1) .

Householder Transformations

Suppose that is a symmetric matrix. Then a sequence of transformations of the form will reduce to a symmetric tridiagonal matrix. Let us visualize the process when . The first transformation is defined to be , where is constructed by applying the above Corollary with the vector being the first column of the matrix . The general form of is

where the letter stands for some element in . As a result, the transformation does not affect the element of :

(2) .

The element denoted is changed because of premultiplication by , and is changed because of postmultiplication by ; since is symmetric, we have . The changes to the elements denoted have been affected by both premultiplication and postmultiplication. Also, since is the first column of , equation (1) implies that .

The second Householder transformation is applied to the matrix defined in (2) and is denoted , where is constructed by applying the Corollary with the vector being the second column of the matrix . The form of is

where stands for some element in . The identity block in the upper-left corner ensures that the partial tridiagonalization achieved in the first step will not be altered by the second transformation . The outcome of this transformation is

(3) .

The elements and were affected by premultiplication and postmultiplication by . Additional changes have been introduced to the other elements by the transformation. The third Householder transformation, , is applied to the matrix defined in (3) where the Corollary is used with being the third column of . The form of is

Again, the identity block ensures that does not affect the elements of , which lie in the upper corner, and we obtain

.

Thus it has taken three transformations to reduce to tridiagonal form.

In general, for efficiency, the transformation is not performed in matrix form. The next result shows that it is more efficiently carried out via some clever vector manipulations.

Theorem (Computation of One Householder Transformation)

If is a Householder matrix, the transformation is accomplished as follows. Let

Let and compute

and
,

then .

Reduction to Tridiagonal Form

Suppose that is a symmetric matrix. Start with

.

Construct the sequence of Householder matrices, so that

for ,

where has zeros below the subdiagonal in columns . Then is a symmetric tridiagonal matrix that is similar to . This process is called Householder's method.

QR method

2010-02-01T08:24:00.000-08:00

Suppose that A is a real symmetric matrix. Householder’s method is used to construct a similar tridiagonal matrix. Then the QR method is used to find all eigenvalues of the tridiagonal matrix. In the latter construction, plane rotations similar to those that were introduced in Jacobi’s method are used to construct the orthogonal matrices . The important step the QR method is the factorization and iteration .

Definition (QR Decomposition): For a nonsingular square matrix , there exists a factorization

where is a unitary and is an upper triangular matrix.

Remark: is a unitary means .

For a real nonsingular square matrix , there exists a factorization

where is an orthogonal matrix and is an upper triangular matrix.

Remark: is a orthogonal means (also ).

QR transformation

After finding the factorization, the transformation is

.

We now investigate a well known and efficient method for finding all the eigenvalues of a general real matrix. The method can be used for an arbitrary real matrix, but for a general matrix it takes many iterations and becomes time consuming.

The method works much faster on special matrices, preferably:
(i) symmetric tri-diagonal,
(ii) Hessenberg matrices,
(iii) symmetric band matrices.

For this module, we will illustrate the QR method for real symmetric matrices.

When solving for eigenvalues of a dense symmetric matrix , the standard practice is to reduce the dense matrix to a tridiagonal matrix through a series of orthogonal transformations, and then to apply an eigenvalue solver for tridiagonal matrices to . The transformations applied to preserve eigenvalues, and the eigenvalues of and are the same.

The popular eigenvalue solver for symmetric tridiagonal matrices is called the implicit method. It applies a series of orthogonal transformations to a tridiagonal matrix which converges to a diagonal matrix . Furthermore, has the same eigenvalues as which are the diagonal elements of . In addition, the product of the orthogonal transformations is a matrix whose columns are the eigenvectors of . The method is called because in each iteration the factorization is computed. The LAPACK routine implementing the implicit algorithm on tridiagonal symmetric matrices is called DSTEQR.

QR Algorithm: The pseudocode for the QR method is:

1. i = 0
2.
3. repeat
4. Factor
5.
6. i = i+1
7. until convergence

Compartment Model

2010-02-01T08:21:00.000-08:00

The Compartment Model I:

The "compartment" model is used to describe the concentration of a dissolved substance in several compartments in a system. For example, the "three-stage system" consists of three tanks containing gallons of solute. Pure water flows at the rate r into the first tank and is mixed, then flows at the rate r into the second tank and is mixed, then flows at the rate r into the third tank and is mixed, and finally flows out of the third tank at the rate r. We define the constants . Then the differential equation for the system is

The Compartment Model II:

Another "compartment" model is used to describe the concentration of a dissolved substance in several compartments in a system. For example, the "closed" three-stage system consists of three tanks containing gallons of solute. The liquid flows at the rate r from the first tank and is mixed, then flows at the rate r into the second tank and is mixed, then flows at the rate r into the third tank and is mixed, and finally flows out of the third tank and back into the first tank at the rate r. We define the constants . Then the differential equation for the system is

The Compartment Model III:

The "five-stage system" consists of five tanks containing gallons of solute. Pure water flows at the rate into the first tank and is mixed, then flows at the rate into the second tank and is mixed, then flows at the rate into successive tanks and is mixed, and finally flows out of the fifth tank at the rate .
We define the constants , , , , and . Then the differential equation for the system is

with the initial conditions

If the size of the tanks is increasing then ordinary eigenvectors are obtained to form the solution. Otherwise, the system might require "generalized eigenvectors" which are important, but we do not want to digress in that direction in this module.

Earthquake Model

2010-02-01T08:16:00.000-08:00

In the study of earthquake induced vibrations on multistory buildings, the free transverse oscillations satisfy the equation

(1) ,

where the forces acting on the i-th floor are

.

Consider a building with n floors each of mass m slugs and the horizontal restoring force of k tons/foot between floors. Then this system reduces to the form

(2) ,
where
.

A horizontal earthquake oscillation of amplitude of the form will produce an acceleration , and the opposite internal force on each floor of the building is . The resulting non-homogeneous system is

, where .

Example 1. Consider a building with n = 6 floors each of mass m = 1250 slugs (weight of 20 tons)
and the horizontal restoring force of k = 10,000 lb/ft = 5 tons/foot between floors.
Then , and this system reduces to the form

Compute the eigenvalues of matrix ,
and the natural frequencies and periods P of oscillation of the building.

Example 2. Solving the above non-homogeneous system for the coefficient vector v for X[t].
The vector v is the solution to the equation .
Use the earthquake amplitude e = 0.075 ft = 0.9 in. for this example.

Solve the linear system using the parameters and e = 0.075.

Find the coefficient vector v and the vector X[t]. Plot the vibrations of each floor.

Matrix Exponential

2010-02-01T08:13:00.000-08:00

We seek a solution of a homogeneous first order linear system of differential equations. For illustration purposes we consider the case:

First, write the system in vector and matrix form

.

Then, find the eigenvalues and eigenvectors of the matrix , denote the eigenpairs of A by

and .

Assumption. Assume that there are two linearly independent eigenvectors , which correspond to the eigenvalues , respectively. Then two linearly independent solution to are

, and
.

Definition (Fundamental Matrix Solution) The fundamental matrix solution , is formed by using the two column vectors .

(1) .

The general solution to is the linear combination

(2) .

It can be written in matrix form using the fundamental matrix solution as follows

.

Notation. When we introduce the notation

,
and

The fundamental matrix solution can be written as

(3) .
or
(4) .

The initial condition

If we desire to have the initial condition , then this produces the equation

.

The vector of constant can be solved as follows

.

The solution with the prescribed initial conditions is

.

Observe that where is the identity matrix. This leads us to make the following important definition

Definition (Matrix Exponential) If is a fundamental matrix solution to , then the matrix exponential is defined to be

.

Notation. This can be written as

(5) ,
or
(6) .

Fact. For a system, the initial condition is

,

and the solution with the initial condition is

,
or
.

Theorem (Matrix Diagonalization) The eigen decomposition of a square matrix A is

,

which exists when A has a full set of eigenpairs for , and d is the diagonal matrix

and

is the augmented matrix whose columns are the eigenvectors of A.

.

Faddeev-Leverrier Method

2010-02-01T08:08:00.000-08:00

Let be an n × n matrix. The determination of eigenvalues and eigenvectors requires the solution of

(1)

where is the eigenvalue corresponding to the eigenvector . The values must satisfy the equation

(2) .

Hence is a root of an nth degree polynomial , which we write in the form

(3) .

The Faddeev-Leverrier algorithm is an efficient method for finding the coefficients of the polynomial . As an additional benefit, the inverse matrix is obtained at no extra computational expense.

Recall that the trace of the matrix , written , is

(4) .

The algorithm generates a sequence of matrices and uses their traces to compute the coefficients of ,

(5)

Then the characteristic polynomial is given by

(6) .

In addition, the inverse matrix is given by

(7) .

Hessenberg Factorization

2010-02-01T08:01:00.000-08:00

An matrix with for is called a Hessenberg matrix. The form of a Hessenberg matrix is

Unitary Matrix

(i) For real matrices, a unitary matrix is a matrix for which .

(ii) For complex matrices, a unitary matrix is a matrix for which .

Hessenberg Factorization of a Matrix:

There are two cases to consider.

(iii) Given a real matrix , there exists a unitary matrix and Hessenberg matrix so that

.

(iv) Given a complex matrix , there exists a unitary matrix and Hessenberg matrix so that

.

Hessenberg Factorization of a Symmetric Matrix:

Given a real symmetric matrix , there exists a unitary matrix and tri-diagonal symmetric matrix so that

.

Remark. This is the case that is easiest to illustrate in a first course in numerical methods.

Golden Ratio

2010-02-01T07:28:00.000-08:00

The function is unimodal on , if there exists a unique number such that

is decreasing on ,
and
is increasing on .

Golden Ratio:

If is known to be unimodal on , then it is possible to replace the interval with a subinterval on which takes on its minimum value. One approach is to select two interior points . This results in . The condition that is unimodal guarantees that the function values and are less than .

If , then the minimum must occur in the subinterval , and we replace b with d and continue the search in the new subinterval . If , then the minimum must occur in the subinterval , and we replace a with c and continue the search in the new subinterval . These choices are shown in Figure 1 below.

If , then squeeze from the right and use the

new interval and the four points .

If , then squeeze from the left and use the

new interval and the four points .

Figure 1. The decision process for the golden ratio search.

The interior points c and d of the original interval , must be constructed so that the resulting intervals , and are symmetrical in . This requires that , and produces the two equations

(1) ,
and
(2) ,

where (to preserve the ordering ).

We want the value of r to remain constant on each subinterval. If r is chosen judicially then only one new point e (shown in green in Figure 1) needs to be constructed for the next iteration. Additionally, one of the old interior points (either c or d) will be used as an interior point of the next subinterval, while the other interior point (d or c) will become an endpoint of the next subinterval in the iteration process. Thus, for each iteration only one new point e will have to be constructed and only one new function evaluation , will have to be made. As a consequence, this means that the value r must be chosen carefully to split the interval of into subintervals which have the following ratios:

(3) and ,
and
(4) and .

If and only one new function evaluation is to be made in the interval , then we must have

.

Use the facts in (3) and (4) to rewrite this equation and then simplify.

,

,

,

.

Now the quadratic equation can be applied and we get

.

The value we seek is and it is often referred to as the "golden ratio." Similarly, if , then it can be shown that .

Fibonacci Method

2010-02-01T07:15:00.001-08:00

An approach for finding the minimum of in a given interval is to evaluate the function many times and search for a local minimum. To reduce the number of function evaluations it is important to have a good strategy for determining where is to be evaluated. Two efficient bracketing methods are the golden ratio and Fibonacci searches. To use either bracketing method for finding the minimum of , a special condition must be met to ensure that there is a proper minimum in the given interval.

The function is unimodal on , if there exists a unique number such that

is decreasing on ,
and
is increasing on .

In the golden ratio search two function evaluations are made at the first iteration and then only one function evaluation is made for each subsequent iteration. The value of remains constant on each subinterval and the search is terminated at the subinterval, provided that or where are the predefined tolerances. The Fibonacci search method differs from the golden ratio method in that the value of is not constant on each subinterval. Additionally, the number of subintervals (iterations) is predetermined and based on the specified tolerances.

The Fibonacci search is based on the sequence of Fibonacci numbers which are defined by the equations

for

Thus the Fibonacci numbers are

Exploration for Fibonacci Numbers

Assume we are given a function that is unimodal on the interval . As in the golden ratio search a value is selected so that both of the interior points will be used in the next subinterval and there will be only one new function evaluation.

If , then the minimum must occur in the subinterval , and we replace and and continue the search in the new subinterval . If , then the minimum must occur in the subinterval , and we replace and and continue the search in the new subinterval . These choices are shown in Figure 1 below.

If , then squeeze from the right and
use the new interval .

If , then squeeze from the left and
use the new interval .

Figure 1. The decision process for the Fibonacci ratio search.

If and only one new function evaluation is to be made in the interval , then we select for the subinterval . We already have relabeled

and since we will relabel it by

then we will have

(1) .

The ratio is chosen so that and and subtraction produces

(2)

And the ratio is chosen so that

(3) .

Now substitute (2) and (3) into (1) and get

(4) .

Also the length of the interval has been shrunk by the factor . Thus and using this in (4) produces

(5) .

Cancel the common factor in (5) and we now have

(6) .

Solving (6) for produces

(7) .

Now we introduce the Fibonacci numbers for the subscript . In equation (7), substitute and get the following

Reasoning inductively, it follows that the Fibonacci search can be begun with

and

for .

Note that the last step will be

,

thus no new points can be added at this stage (i.e. the algorithm terminates). Therefore, the set of possible ratios is

.

There will be exactly n-2 steps in a Fibonacci search!

The subinterval is obtained by reducing the length of the subinterval by a factor of . After steps the length of the last subinterval will be

.

If the abscissa of the minimum is to be found with a tolerance of, then we want . It is necessary to use n iterations, where n is the smallest integer such that

(8) .

Note. Solving the above inequality requires either a trial and error look at the sequence of Fibonacci numbers, or the deeper fact that the Fibonacci numbers can be generated by the formula

.

Knowing this fact may be useful, but we still need to compute all the Fibonacci numbers in order to calculate the ratios .

The interior points and of the subinterval are found, as needed, using the formulas

(9) ,

(10) .

Note. The value of n used in formulas (9) and (10) is found using inequality (8).

Quadratic and Cubic Methods

2010-02-01T07:03:00.000-08:00

The function is unimodal on , if there exists a unique number such that

is decreasing on ,
and
is increasing on .

Minimization Using Derivatives

Suppose that is unimodal over and has a unique minimum at . Also, assume that is defined at all points in . Let the starting value lie in . If , then the minimum point p lies to the right of . If , then the minimum point p lies to the left of .

Our first task is to obtain three test values,

(1) ,
so that
(2) .

Suppose that ; then and the step size h should be chosen positive. It is an easy task to find a value of h so that the three points in (1) satisfy (2). Start with in formula (1) (provided that ); if not, take , and so on.

Case (i) If (2) is satisfied we are done.

Case (ii) If , then .
We need to check points that lie farther to the right. Double the step size and repeat the process.

Case (iii) If , we have jumped over p and h is too large.
We need to check values closer to . Reduce the step size by a factor of and repeat the process.

When , the step size h should be chosen negative and then cases similar to (i), (ii), and (iii) can be used.

Quadratic Approximation to Find p

Finally, we have three points (1) that satisfy (2). We will use quadratic interpolation to find , which is an approximation to p. The Lagrange polynomial based on the nodes in (1) is

(3),

where .

The derivative of is

(4).

Solving in the form yields

(5).

Multiply each term in (5) by and collect terms involving :

This last quantity is easily solved for :

.

The value is a better approximation to p than . Hence we can replace with and repeat the two processes outlined above to determine a new h and a new . Continue the iteration until the desired accuracy is achieved. In this algorithm the derivative of the objective function was used implicitly in (4) to locate the minimum of the interpolatory quadratic. The reader should note that the subroutine makes no explicit use of the derivative.

Cubic Approximation to Find p

We now consider an approach that utilizes functional evaluations of both and . An alternative approach that uses both functional and derivative evaluations explicitly is to find the minimum of a third-degree polynomial that interpolates the objective function at two points. Assume that is unimodal and differentiable on , and has a unique minimum at . Let . Any good step size h can be used to start the iteration. The Mean Value Theorem could be used to obtain and if was just to the right of the minimum, then the slope might be twice which would mean that we do not know how much further to the right lies, so we can imagine that is close to and estimate h with the formula:

.

Thus . The cubic approximating polynomial is expanded in a Taylor series about (which is the abscissa of the minimum). At the minimum we have , and we write in the form:

(6),
and
(7).

The introduction of in the denominators of (6) and (7) will make further calculations less tiresome. It is required that , , , and . To find we define:

(8) ,

and we must go through several intermediate calculations before we end up with .

Use use (6) to obtain

Then use (8) to get

Then substitute and we have

(9)

Use use (7) to obtain

Then use (8) to get

Then substitute and we have

(10)

Finally, use (7) and write

Then use (8) to get

(11)

Now we will use the three nonlinear equations (9), 10), (11) listed below in (12). The order of determining the variables will be (the variable will be eliminated).

(12)

First, we will find which is accomplished by combining the equation in (12) as follows:

Straightforward simplification yields , therefore is given by

(13) .

Second, we will eliminate by combining the equation in (12) as follows, multiply the first equation by and add it to the third equation

which can be rearranged in the form

Now the quadratic equation can be used to solve for

It will take a bit of effort to simplify this equation into its computationally preferred form.

Hence,

(14)

Therefore, the value of is found by substituting the calculated value of in (14) into the formula . To continue the iteration process, let and replace and with and , respectively, in formulas (12), (13), and (14). The algorithm outlined above is not a bracketing method. Thus determining stopping criteria becomes more problematic. One technique would be to require that , since .