TI83 Calculator Programs for Numerical Analysis Problems - Part 1

These programs are copyrighted (1997-2007), but you may copy them for instructional purposes as long as no profit is made from their use. The author is not responsible for any data loss which may be caused to any calculator or its memory by the use of these programs.

The following are TI-83 calculator programs for the solution of numerical analysis problems of a type which are usually found in college courses or textbooks in numerical analysis. These programs are also suitable, with minor modifications, for use in other TI calculators, such as the TI82, TI85, TI86, TI89, TI92, etc... . (The easier programs could be used, with minor modifications, in a TI81 calculator.) These TI-83 programs, together with their descriptions, will probably be read more easily by students who are either taking a course in numerical analysis, or who have already taken a course in numerical analysis. (Click here for comments concerning calculators other than the TI83.)

Comments and questions concerning these programs may be addressed to Gerald Roskes, Department of Mathematics, Queens College, Flushing, New York 11367, or send email to gerald_roskes@qc.edu.

This web page will be continually expanding as we add more programs to our list. If you have an interest in calculator programs for numerical analysis, you should view our page every few months. (This web page was last updated in December, 2007.)

The following are links to the various TI83 programs:

Part 1 (Programs 1 - 23)

Prog1 SIMPITER (Simple Iteration)
Prog2 BISECT (Bisection Method)
Prog3 SECANT (Secant Method)
Prog4 STEFFEN (Steffensen's Method)
Prog5 AITKEN (Aitken's Delta2 Method)
Prog6 HORNER (Horner's Method)
Prog7 FALSEPOS (False Position Method)
Prog8 MULLER (Muller's Method)
Prog9 SYSITER2 (Iteration for 2 x 2 Systems)
Prog10 SYSNEWT2 (Newton Iteration for 2 x 2 Systems)
Prog11 JACOBI (Jacobi Iteration)
Prog12 GSEIDEL (Gauss-Seidel Iteration)
Prog13 SOR (Successive Over-Relaxation)
Prog14 SCLPWR (Scaled Power Method)
Prog15 WDEFLATE (Wielandt Deflation)
Prog16 INTERPN (N-Point Interpolation)
Prog17 INTERP2 (2-Point Interpolation)
Prog18 INTERP3 (3-Point Interpolation)
Prog19 INTERP4 (4-Point Interpolation)
Prog20 INTERP5 (5-Point Interpolation)
Prog21 NEVILLE (Neville's Method)
Prog22 NEWTDIV (Newton's Divided Difference Method)
Prog23 LAGRANGE (Lagrange Interpolation Polynomial)

Part 2 (Programs 24 - 42)

Prog24 HERMITEN (N Point Hermite Interpolation)
Prog25 HERMITDD (Hermite Divided Difference Method)
Prog26 HERMITLA (Hermite-Lagrange Interpolating Polynomial)
Prog27 DERIV2PT (Derivative 2-Point Approximation)
Prog28 DERIV3PT (Derivative 3-Point Approximation)
Prog29 DERIV5PT (Derivative 5-Point Approximation)
Prog30 DDER3PT (Double Derivative 3-Point Approximation)
Prog31 NEWFORDF (Newton's Forward Difference Method)
Prog32 INEWCOTE (Newton-Cotes Integration Method
Prog33 ICOMTRAP (Composite Trapezoid Rule)
Prog34 ICOMMIDP (Composite Midpoint Rule)
Prog35 ICOMSIMP (Composite Simpson Rule)
Prog36 IROMBERG (Romberg Integration)
Prog38 IGAUSSAB (Gaussian Quadrature Method on [a,b])
Prog39 IGAUSAB2 (Gauss 2 Point Quadrature on [a,b])
Prog40 IGAUSAB3 (Gauss 3 Point Quadrature on [a,b])
Prog41 IGAUSAB4 (Gauss 4 Point Quadrature on [a,b])
Prog42 IGAUSAB5 (Gauss 5 Point Quadrature on [a,b])

Part 3 (Programs 43 - 52)

Prog43 DEEULER (Euler's Method)
Prog44 DETAYLR2 (Taylor Method, Order 2)
Prog45 DETAYLR4 (Taylor Method, Order 4)
Prog46 DEMIDPT (DE Midpoint Method)
Prog47 DEMODEUL (DE Modified Euler Method)
Prog48 DEHEUN (DE Heun's Method)
Prog49 DERK4 (DE Runge Kutta 4th Order Method)
Prog50 DEBASH2 (Adams Bashforth 2 Step)
Prog51 DEMOULT2 (Adams Moulton 2 Step)
Prog52 DEPCB2M2 (Bashforth 2 Step Predictor, Moulton 2 Step Corrector)

Part 4 (Programs 53 - 70)

Prog53 DEBASH3 (Adams Bashforth 3 Step)
Prog54 DEMOULT3 (Adams Moulton 3 Step)
Prog55 DEPCB4M3 (Bashforth 4 Step Predictor, Moulton 3 Step Corrector)
Prog56 DERKV6 (DE Runge Kutta Verner 6th Order Method)
Prog57 LSPOLY1D (Least Squares Linear Polynomial (Discrete Data))
Prog58 LSPOLY2D (Least Squares Quadratic Polynomial (Discrete Data))
Prog59 LSPOLY3D (Least Squares Cubic Polynomial (Discrete Data))
Prog60 LSPOLY4D (Least Squares Quartic Polynomial (Discrete Data))
Prog61 LSEXPD (Least Squares Exponential Function (Discrete Data))
Prog62 LSPOLY1C (Least Squares Linear Polynomial (Continuous Data))
Prog63 LSPOLY2C (Least Squares Quadratic Polynomial (Continuous Data))
Prog64 LSPOLY3C (Least Squares Cubic Polynomial (Continuous Data))
Prog65 LSPOLY4C (Least Squares Quartic Polynomial (Continuous Data))
Prog66 LSLEGNDR (Least Squares Legendre Polynomials ( Degrees 0-5, [ -1 , 1 ] ) )
Prog67 LSLGDRAB (Least Squares Legendre Polynomials ( Degrees 0-5, [ a , b ] ) )
Prog68 SYSNEWT3 (Newton Iteration for 3 x 3 Systems)
Prog69 ECONOMIZ (Chebyshev Economization, (Degrees 0-8, [ -1 , 1 ] ) )
Prog70 ECONOMAB (Chebyshev Economization, (Degrees 0-8, [ a , b ] ) )

Part 5 (Programs 71 - 86)

Prog71 OPTINTP2 (Optimal Interolation, 2 Points)
Prog72 OPTINTP3 (Optimal Interolation, 3 Points)
Prog73 OPTINTP4 (Optimal Interolation, 4 Points)
Prog74 OPTINTP5 (Optimal Interolation, 5 Points)
Prog75 LSCHEBYS (Least Squares Chebyshev Polynomials, (Degrees 0-3, [ -1 , 1 ] ) )
Prog76 LSCHEBAB (Least Squares Chebyshev Polynomials, (Degrees 0-3, [ a , b ] ) )
Prog77 NEWTON (Newton Iteration with Test)
Prog78 MODNEWT (Modified Newton Method)
Prog79 UNSCLPWR (Unscaled Power Method)
Prog80 BEZIER (Bezier Curve, 2 Points)
Prog81 LSTRCN2 (LS Trig Curve , Continuous Data , Case N=2)
Prog82 LSTRCN3 (LS Trig Curve , Continuous Data , Case N=3)
Prog83 LSTRDN2 (LS Trig Curve , Discrete Data , Case N=2 , M > 2)
Prog84 LSTRDN3 (LS Trig Curve , Discrete Data , Case N=3 , M > 3)
Prog85 LSTRIN2 (LS Trig Interpolation Curve , Case N=2)
Prog86 LSTRIN2 (LS Trig Interpolation Curve , Case N=3)
Prog87 SYSITER3 (Iteration for 3 x 3 Systems)
Prog88 SYSGS2 (Gauss-Seidel for 2 x 2 Nonlinear Systems)
Prog89 SYSGS3 (Gauss-Seidel for 3 x 3 Nonlinear Systems)

-------------------------------------------------------------------------

Program #1 SIMPITER (Simple Iteration)

This is our shortest and simplest program. It has only one instruction. The iteration is Xn+1 = f ( Xn ) . Before executing the program, the function f(X) is stored in the function variable Y1 , and the initial value X0 is stored in the variable X. Next, execute the program. After the program is executed, keep pressing ENTER. The iterates will appear on the home screen.

In order to enter the term "Y1" into the program, press [VARS], [right arrow], [ENTER], [ENTER]. Notice that we do not use the "Stop" or "Return" instruction as a second instruction. These instructions would not display any result on the home screen. Also, we do not put any blank lines into the program after instruction#1, as this would cause the calculator to display only the message "Done" on the home screen, and nothing more.

Instruction #1_____Y1 -> X..........(Note: We use the symbol "->" to
...............................................................represent the operation "Store".)

The following is a program that displays 5 iterates with each push of the ENTER button. Before executing the program, put f(X) into Y1, and put X0 into X. Then execute the program. Keep pressing ENTER.

Instruction#1_____Disp "X ITERATES ARE"
Instruction#2_____For(I,1,4)
Instruction#3_____Y1 -> X
Instruction#4_____Disp X
Instruction#5_____End
Instruction#6_____Y1 -> X

-------------------------------------------------------------------------

Program #2 BISECT (Bisection Method)

Let f(X) be defined on [a,b]. We require that f(X) be continuous, and that the product f(a)*f(b) is negative. Then f(X) has at least one zero in the interval [a,b]. The bisection method will converge to one of these zeroes. This program will compute and display the iterates in the bisection method. Before executing the program, the numbers a and b are stored in the variables A and B respectively, and the function f(X) is stored in Y1. Next, execute the program. The program computes the midpoint X = (A+B)/2 (instruction #3) and displays A,B,X, f(X). If f(B)*f(X) is positive (so that f(B) and f(X) are of the same sign), then B is replaced by the midpoint X (instruction 8). Otherwise, A is replaced by the midpoint. After executing the program, keep pressing ENTER. The iterates will appear on the home screen.

Instruction #1_____B -> X..........(Note: We use the symbol "->" to
Instruction #2_____Y1 -> D.........represent the operation "Store".)
Instruction #3_____(A+B)/2 -> X
Instruction #4_____Y1 -> Q
Instruction #5_____Disp A,B,X,Q
Instruction #6_____If D*Q > 0
Instruction #7_____Then
Instruction #8_____X -> B
Instruction #9_____Else
Instruction #10____X -> A
Instruction #11____End

-------------------------------------------------------------------------

Program #3 SECANT (Secant Method)

In the secant method, we start with a function f(X), and two values X0, X1 close to a zero of f(X). (The values X0 and X1 should be different.) The points (X0,f(X0)) and (X1,f(X1))determine a line. In general, this line will intersect the X-axis. The point of intersection with the X-axis determines X2. The values X1,X2 are then used to update the values X0, X1, and the process is repeated.

Before executing the program, we store X0 in the variable A, X1 in the variable B, and f(X) in the function variable Y1. Next, execute the program. When the program is executed, the values of f(X0) and f(X1) are stored in the variables C and D (see instructions #2 and #4). The value of X2 is computed from X2 = X1 - f(X1)*(X1 - X0) / (f(X1) - f(X0)) (see instruction #6).

In general, we get rapid convergence to a zero of f(X). In this case, the program will eventually encounter a zero denominator in instruction #6 when, to machine accuracy, f(X1) = f(X2) = 0. We test for the condition that f(X1) = f(X2) in instruction#5. If this condition is satisfied, then X2 is left equal to X1 (from instruction #3). (In this case, one should further check that f(X1) = 0 to be sure that we have convergence to a zero of f(X).)

When the program is executed, X2 is displayed (from instruction #8), and X0, X1 are updated (instructions #7,8). After execution, keep pressing ENTER. The iterates appear on the home screen.

(Our last instruction #8 produces a result (namely, X2) which is also displayed on the home screen. Notice that we do not use the "Stop" or "Return" instruction as our last instruction. These instructions would not display any result on the home screen.)

Instruction #1_____A -> X..........(Note: We use the symbol "->" to
Instruction #2_____Y1 -> C.........represent the operation "Store".)
Instruction #3_____B -> X
Instruction #4_____Y1 -> D
Instruction #5_____If D - C != 0____________(See note below.)
Instruction #6_____B - D(B - A) / (D - C) -> X
Instruction #7_____B -> A
Instruction #8_____X -> B

Note that in instruction#5, we are using the symbol != to represent "not equal". When entering this program into the TI83, one should use the more standard symbol from the TEST menu.

-------------------------------------------------------------------------

Program #4 STEFFEN (Steffensen's Method)

We assume we have an iteration equation X = g(X) and a value X0 near a fixed point of g(X). If we set A = X0, B = g(A), C = g(B), then Steffensen's method computes X1 from Aitken's equation: X1 = A - (B -A)2 / (C - 2 B + A). The number X1 is used to update X0, and the process is repeated. In general, Steffensen's method will converge quadratically to a solution of X = g(X), even though the equation X = g(X) produces linearly convergent, or even divergent, iterates.

Before executing the program, we store X0 in the variable X, and the function g(X) in the function variable Y1. The program is then executed, and will display A, B, C, and X1 (see instruction#6). The program also updates X0 with X1. Because of convergence, we will eventually obtain values for A, B, C which are all equal to X0 (to machine accuracy). This will create a zero denominator in Aitken's equation. If the denominator in Aitken's equation vanishes (see instruction#4), then X1 is left equal to X0.

After executing the program, keep pressing ENTER. The Steffensen iterates (for A, B, C, X1) will appear on the home screen.

Instruction #1_____X -> A..........(Note: We use the symbol "->" to
Instruction #2_____Y1 -> B.........represent the operation "Store".)
Instruction #3_____Y1(B) -> C
Instruction #4_____If C - 2 B + A != 0____________(See note below.)
Instruction #5_____X - ( B - A )2 / ( C - 2 B + A ) -> X
Instruction #6_____Disp A,B,C,X

Note that in instruction#4, we are using the symbol != to represent "not equal". When entering this program into the TI83, one should use the more standard symbol from the TEST menu.

-------------------------------------------------------------------------

Program #5 AITKEN (Aitken's Delta2 Method)

We assume {Xn} is generated from the iteration Xn+1 = g(Xn). We also assume Xn converges linearly to a number S. Aitken's method produces a new sequence {Zn} which converges more quickly to S. We define A = Xn, B = Xn+1, C = Xn+2 . Then the Aitken equation Zn = A - (B - A)2 / (C - 2B + A) defines the new sequence. Our program displays Xn and Zn so that the more rapid comvergence of Zn can be observed.

Before executing our program, we store X0 in the variable X, and g(X) in the function variable Y1. The program is then executed. The value of Z0 is computed from Aitken's equation (instruction#5). The program displays X0 and Z0 (instruction#7), and updates X0 with X1 (instruction#6).

After executing the program, keep pressing ENTER. The X and Z iterates will be displayed on the home screen.

Since Xn converges to S, the values of A, B, C will eventually all equal to S to machine accuracy. This creates a zero denominator in instruction#5. We test for this condition (instruction#4). When this condition occurs, the value of Zn is left equal to Zn-1 , which, in general, will be S to machine accuracy.

Instruction #1_____X -> A..........(Note: We use the symbol "->" to
Instruction #2_____Y1 -> B.........represent the operation "Store".)
Instruction #3_____Y1(B) -> C
Instruction #4_____If C - 2 B + A != 0____________(See note below.)
Instruction #5_____X - ( B - A )2 / ( C - 2 B + A ) -> Z
Instruction #6_____B -> X
Instruction #7_____Disp A,Z

Note that in instruction#4, we are using the symbol != to represent "not equal". When entering this program into the TI83, one should use the more standard symbol from the TEST menu.

-------------------------------------------------------------------------

Program #6 HORNER (Horner's Method)

In Horner's method, we assume we have a polynomial of degree n given by P(X) = anXn + an-1 Xn-1 + ... + a0 . We also assume we are given a number X0, possibly close to a zero of P(X). Horner's method computes two sequences as follows: bn = an, bi = bi+1 X0 + ai , i = n-1, n-2, ... , 0 , and cn = bn , ci = ci+1 X0 + bi , i = n-1, n-2, ... , 0 . Then it can be shown that P(X0) = b0 and P'(X0) = c1. (Although we compute c0, it is not needed for our purposes.) Our program computes the b,c sequences and the Newton iterate X1 = X0 - P(X0) / P'(X0) . The value of X0 is then updated by X1.

Before executing our program, we first put X0 into the variable X, and the list { an, an-1, ... , a0 } into the list variable L1. (Thus, L1(1) = an .) The program is then executed. Lists L2, L3 are dimensioned to store n+1 numbers (instructions#2,3). Horner's method is then used to compute (see instructions#4,5,6,7,8,9) the list {bn, bn-1, ... , b0}, (which is stored in L2) and the list {cn, cn-1, ... , c0} (which is stored in L3). The Newton iterate X1 is computed and used to update X0 (instruction#10). X1 is also put into L3(N+1), replacing the unneeded c0. For output purposes, L2 and L3 are put into the columns of a matrix [A] (instruction#12), and the matrix is displayed (instruction#13). Then, P(X0) = [A](N+1,1) , P'(X0) = [A](N,2) , and X1 = [A](N+1,2) . Keep pressing ENTER. Convergence is rapid if X0 is close to a zero of P(X).

Instruction #1_____dim(L1) - 1 -> N..........(Note: We use the symbol "->" to
Instruction #2_____N + 1 -> dim(L2).........represent the operation "Store".)
Instruction #3_____N + 1 -> dim(L3)
Instruction #4_____L1(1) -> L2(1)
Instruction #5_____L2(1) -> L3(1)
Instruction #6_____For(I,1,N)
Instruction #7_____X L2(I) + L1(I+1) -> L2(I+1)
Instruction #8_____X L3(I) + L2(I+1) -> L3(I+1)
Instruction #9_____End
Instruction #10____X - L2(N+1) / L3(N) -> X
Instruction #11____X -> L3(N+1)
Instruction #12____List>matr(L2 , L3 , [A] ).........(See note below.)
Instruction #13____[A]

Note: Use the MATRIX menu to enter the matrix name [A] onto a program line. Do not type the characters "[" , "A" , "]" separately.

-------------------------------------------------------------------------

Program #7 FALSEPOS (False Position Method)

Let f(X) be defined on [a,b]. We require that f(X) be continuous, and that the product f(a)*f(b) is negative. Then f(X) has at least one zero in the interval [a,b]. The false position method will converge to one of these zeroes. This program will compute and display the iterates in the false position method.

Before executing the program, the numbers a and b are stored in the variables A and B respectively, and the function f(X) is stored in Y1 . Next, execute the program. The program computes the zero of the secant line joining the points (A,f(A)) and (B,f(B)) (see instructions#1,2,3,4,5). This zero is put into the variable X. Thus, X = B - f(B) ( B - A ) / ( f(B) - f(A) ) . The values of A,B,X, and f(X) are displayed (see instruction#7). If f(B) * f(X) is positive (so that f(B) and f(X) are of the same sign), then B is replaced by X (see instruction#10). Otherwise, A is replaced by X (see instruction#12). After executing the program, keep presing ENTER. The iterates will appear on the home screen. In general, convergence is linear. Note that in this method, the updated values for A and B will, in general, form a new bracket on a zero of f(X).

Instruction#1_____A -> X..........(Note: We use the symbol "->" to
Instruction#2_____Y1 -> C.........represent the operation "Store".)
Instruction#3_____B -> X
Instruction#4_____Y1 -> D
Instruction#5_____B - D ( B - A ) / ( D - C ) -> X
Instruction#6_____Y1 -> Q
Instruction#7_____Disp A,B,X,Q
Instruction#8_____If D * Q > 0
Instruction#9_____Then
Instruction#10____X -> B
Instruction#11____Else
Instruction#12____X -> A
Instruction#13____End

-------------------------------------------------------------------------

Program #8 MULLER (Muller's Method)

Muller's method is used to determine a zero of a polynomial f(X). The iteration starts with three distinct approximations X0,X1,X2 to a zero of f(X). The three data points (X0,f(X0)), (X1,f(X1)), (X2,f(X2)) are used to determine a parabola P(X). The zero X* of the parabola which is closest to X2 is used as the next approximation to the zero of f(X). The values of X0,X1,X2 are updated with the values of X1,X2,X* respectively, and the process is repeated.

In order to reduce roundoff error, the parabola P(X) is written as P(X) = a(X-X2)^2 + b(X-X2) + c . (That is, we "center" the computation at X2.) The three data points determine three equations for a,b,c :
a(X0-X2)^2 + b(X0-X2) + c = f(X0) ,
a(X1-X2)^2 + b(X1-X2) + c = f(X1) ,
c = f(X2) .
Since the TI83 does not support matrices with complex entries, we solve for a,b,c using Cramer's rule. Thus,
a = [(X1-X2)(f(X0)-f(X2)) - (X0-X2)(f(X1)-f(X2))] / [(X0-X2)(X1-X2)(X0-X1)] ,
b = [((X0-X2)^2)(f(X1)-f(X2)) - ((X1-X2)^2)(f(X0)-f(X2))] / [(X0-X2)(X1-X2)(X0-X1)] ,
c = f(X2) .
The zeroes of P(X) are given by the quadratic formula as
X2 + (-b + (b^2 - 4ac)^.5) / (2a) , and
X2 + (-b - (b^2 - 4ac)^.5) / (2a) .
The one closest to X2 is used as X*.

Before executing the program, the valuse of X0,X1,X2 are placed into the variables D,E,F respectively. The polynomial f(X) is placed into Y1. The program computes f(X0),f(X1),f(X2) and stores these values into the variables G,H,I respectively (see instructions#1,2,3). Instruction#4 replaces (temporarily) the contents of D,E with (X0-X2) and (X1-X2) in order to prepare for Cramer's rule. The values of a,b,c are stored into A,B,C respectively (see instructions#14,15,16). The value of X* is computed (see instructions#17,18,28,29,30,31,32) and placed into the variable X. Finally the variables D,E are restored to their original values (instruction#34), and an output matrix [A] is set up (instructions#35,36,37). The values for X0,X1,X2 are updated (instruction#38), and the matrix [A] is displayed (instruction#39). In the output matrix [A], columns 1 and 2 contain the real and imaginary parts of X0,X1,X2,X*. Columns 3 and 4 contain the real and imaginary parts of the corresponding function values from the polynomial f(X).

The program checks for some unlikely events:
1. If either X0,X1 or X2 is a zero of f(X), a message is displayed (see instructions#5,6,7,8). The data is also displayed.
2. If X0=X1 or X1=X2 or X0=X2, a message is displayed (see instructions# 9,10,11,12). The data is also displayed. In this case, new and distinct values for X0,X1,X2 should be chosen, and the program restarted.
3. If the data points are on a horizontal line, a message is displayed (see instructions#19,20,21,22,23). In this case, the program stops, and new values for X0,X1,X2 should be chosen.
4. If the data is on an oblique line, a message is displayed, and the program continues (see instructions#24,25,26,27). X* is computed as the intersection point of the line with the X-axis.

As mentioned above, the values of X0,X1,X2 are placed into the variables D,E,F before executing the program. Also, MODE should be set to complex (this setting is "a+bi" in the MODE menu.) After executing the program, keep pressing ENTER. The iterates will appear on the home screen. (You may need to scroll through the output matrix [A] on the home screen.) In general, convergence is rapid.

(See also the short version of Muller's program toward the end of this section.)

Instruction#1_____D -> X : Y1 -> G..........(Note: We use the symbol "->" to
Instruction#2_____E -> X : Y1 -> H...........represent the operation "Store".)
Instruction#3_____F -> X : Y1 -> I
Instruction#4_____D - F -> D : E - F -> E
Instruction#5_____If abs(GHI) = 0
Instruction#6_____Then : Disp "ZERO APPEARS IN" : Disp "FOLLOWING TABLE"
Instruction#7_____Go to 1
Instruction#8_____End
Instruction#9_____( D - E ) D E -> J
Instruction#10____If abs ( J ) = 0 : Then
Instruction#11____Disp "DOUBLE X VALUES"
Instruction#12____Go to 1
Instruction#13____End
Instruction#14____(D(I-H)+E(G-I)) / J -> A
Instruction#15____(D2(H-I)+E2(I-G)) / J -> B
Instruction#16____I -> C
Instruction#17____(B2 - 4AC) ^ (1/2) -> S
Instruction#18____2A -> T
Instruction#19____If abs(A)=0 and abs(B)=0
Instruction#20____Then
Instruction#21____Disp "HORIZ LINE" : Disp "CHANGE DATA"
Instruction#22____D + F -> D : E + F -> E
Instruction#23____Stop : End
Instruction#24____If abs(A)=0 and abs(B) !=0_________(See note below)
Instruction#25____Then : F - C / B -> X
Instruction#26____Disp "OBLIQUE LINE"
Instruction#27____Go to 1 : End
Instruction#28____( (-) B + S ) / T -> K
Instruction#29____( (-) B - S ) / T -> L
Instruction#30____If abs(K) <= abs(L)_________(See note below)
Instruction#31____Then : F + K -> X
Instruction#32____Else : F + L -> X : End
Instruction#33____Lbl 1
Instruction#34____D + F -> D : E + F -> E
Instruction#35____{ D, E, F, X } -> L1
Instruction#36____{ G, H, I, Y1 } -> L2
Instruction#37____List>matr ( real (L1) , imag (L1) , real (L2) , imag (L2) , [A] )
Instruction#38____E -> D : F -> E : X -> F
Instruction#39____[A].................................................(See note below.)

Note that in instruction#24, we are using the symbol != to represent "not equal". When entering this program into the TI83, one should use the more standard symbol from the TEST menu. Also, in instruction#30, we are using the symbol <= to represent "less than or equal to". When entering this program into the TI83, one should use the more standard symbol from the TEST menu.
For instructions#37,39, use the MATRIX menu to enter the matrix name [A] onto the program lines. Do not type the characters "[" , "A" , "]" separately.

Short version of Muller's program

For most textbook problems, we can shorten our program and still obtain adequate results. We delete instructions#5,6,7,8 (test for a zero). We also delete instructions#10,11,12,13 (test for duplicate data values). We also delete instructions#19,20,21,22,23 (test for a horizontal line). We also delete instructions#24,25,26,27 (test for an oblique line). Finally, we shorten the output (instructions#33,34,35,36,37,38,39 are replaced) to display just the real and imaginary parts of X* and f(X*). The resulting program is listed below.

After executing the short version, keep pressing ENTER. The iterates will appear on the home screen. The program will stop when it encounters a zero denominator. This usually happens after the iterates have converged to machine accuracy.

Instruction#1_____D -> X : Y1 -> G
Instruction#2_____E -> X : Y1 -> H
Instruction#3_____F -> X : Y1 -> I
Instruction#4_____D - F -> D : E - F -> E
Instruction#5_____( D - E ) D E -> J
Instruction#6_____(D(I-H)+E(G-I)) / J -> A
Instruction#7_____(D2(H-I)+E2(I-G)) / J -> B
Instruction#8_____I -> C
Instruction#9_____(B2 - 4AC) ^ (1/2) -> S
Instruction#10____2A -> T
Instruction#11____( (-) B + S ) / T -> K
Instruction#12____( (-) B - S ) / T -> L
Instruction#13____If abs(K) <= abs(L)_________(See note below)
Instruction#14____Then : F + K -> X
Instruction#15____Else : F + L -> X : End
Instruction#16____E + F -> D : F -> E : X -> F : Y1 -> I
Instruction#17____List>matr( {real (F) , imag (F) , real (I) , imag (I)} , [A] )
Instruction#18____[A]........................(See note below.)

Note that in instruction#13, we are using the symbol <= to represent "less than or equal to". When entering this program into the TI83, one should use the more standard symbol from the TEST menu.
For instructions#17,18, use the MATRIX menu to enter the matrix name [A] onto the program lines. Do not type the characters "[" , "A" , "]" separately.

-------------------------------------------------------------------------

Program #9 SYSITER2 (Iteration for 2 x 2 Systems)

This program solves the system X = f1(X,Z) , Z = f2(X,Z) by iteration. If (X0,Z0) is an approximate solution, we compute (X1,Z1) from the equations X1 = f1(X0,Z0) , Z1 = f2(X0,Z0) . The values for X1 , Z1 are used to update X0 , Z0 , and the process is repeated.

Before executing the program, the function f1(X,Z) is stored in the function variable Y1 (as an expression in X and Z), and the function f2(X,Z) is stored in the function variable Y2 (as an expression in X and Z).The value for X0 is stored in the variable X, and the value for Z0 is stored in the variable Z. Instruction#1 computes X1 and stores the value temporarily in the variable W. Instruction#2 computes Z1 and stores the value in the variable Z (thus,Z is "updated"). Instruction#3 updates X with the value in W (which is X1). Instruction#4 displays the updated values for (X,Z).

After executing the program, keep pressing ENTER. The iterates will appear on the home screen.

Instruction#1_____Y1 -> W..........(Note: We use the symbol "->" to
Instruction#2_____Y2 -> Z............represent the operation "Store".)
Instruction#3_____W -> X
Instruction#4_____Disp X,Z

-------------------------------------------------------------------------

Program #10 SYSNEWT2 (Newton Iteration for 2 x 2 Systems)

This program solves the system f1(X,Z) = 0 , f2(X,Z) = 0 using Newton iteration for systems. If (X0,Z0) is an approximation to a solution, the Newton iteration computes a new approximation (X1 , Z1) from the vector equation X1 = X0 - [J]-1F(X0 , Z0), where the vector X0 = (X0 , Z0)T , the vector X1 = (X1 , Z1)T , the vector function F = (f1 , f2)T and where the matrix [J] is the Jacobian matrix for F at (X0 , Z0) .The values for X1 , Z1 are used to update the values of X0 , Z0 , and the process is repeated.

Before executing the program, we store f1(X , Z) in the function variable Y1 (as an expression in X and Z), and we store f2(X , Z) in the function variable Y2 (as an expression in X and Z). Also, X0 is stored into the variable X , and Z0 is stored into the variable Z. The program computes the components of the vector F(X0 , Z0) in instructions#3,4 and stores the result in the column matrix [B]. (The dimensions of the matrices [A] and [B] are set by instructions#1,2.) The partial derivatives Dxf1 , Dzf1 , Dxf2 , Dzf2 are computed at (X0 , Z0) and stored in the matrix [A] (see instructions#5,6,7,8). Matrix [A] is the Jacobian matrix for the vector function F. The values for X1 , Z1 are computed in instructions#9,10,11 and are used to update X0 , Z0. An output matrix [D] is constructed in instructions#12,13,14. Column 1 of [D] has the vector (X1 , Z1)T. Columns 2,3 of [D] contain the Jacobian matrix [J] at (X0 , Z0). Columns 4,5 of [D] contain the inverse matrix [J]-1 at (X0 , Z0). Instruction#15 displays the output matrix [D].

After executing the program, keep pressing ENTER. The iterates will appear on the home screen. In general, convergence is quadratic to machine accuracy.

Instruction#1_____{2,2} -> dim([A])..........(Note: We use the symbol "->" to
Instruction#2_____{2,1} -> dim([B])...........represent the operation "Store".)
Instruction#3_____Y1 -> [B](1,1)...............(See note below.)
Instruction#4_____Y2 -> [B](2,1)
Instruction#5_____Nderiv(Y1, X, X, .00001) -> [A](1,1)
Instruction#6_____Nderiv(Y1, Z, Z, .00001) -> [A](1,2)
Instruction#7_____Nderiv(Y2, X, X, .00001) -> [A](2,1)
Instruction#8_____Nderiv(Y2, Z, Z, .00001) -> [A](2,2)
Instruction#9_____[A]-1 [B] -> [C]
Instruction#10____X - [C](1,1) -> X
Instruction#11____Z - [C](2,1) -> Z
Instruction#12____Matr>list([A], L1, L2)
Instruction#13____Matr>list([A]-1, L3, L4)
Instruction#14____List>matr( {X , Z} , L1 , L2 , L3 , L4 , [D] )
Instruction#15____[D]

Note: Use the MATRIX menu to enter the matrix names [A] ,[B], etc... onto a program line. Do not type the characters "[" , "A" , "]" , "B" , ... separately.

(Note: To obtain a display of output which includes the
vectors F(X0 , Z0) and [J]-1F(X0 , Z0), replace Instruction#15 with
the instruction: augment ( augment ( [D] , [B] ) , [C] ) -> [D] . )

-------------------------------------------------------------------------

Program #11 JACOBI (Jacobi Iteration)

Jacobi iteration is a method for solving a linear system of equations Ax = b , where A is an n by n matrix. In matrix form, Jacobi iteration may be described as follows:

We "split" the matrix A into D - L - U, where D is a diagonal matrix with Dii = Aii, L is a strictly lower triangular matrix with Lij = -Aij when i > j, and U is a strictly upper triangular matrix with Uij = -Aij when j > i. We require that Aii is not zero for each i. Then, D-1 exists, and AX = b may be rewritten as x = D-1(L + U) x + D-1b. If we set Tj = D-1( L + U ) and cj = D-1 b , then the original system is equivalent to the system x = Tjx + cj . If the spectral radius of Tj is less than 1, then we can solve by iteration. Thus, we choose the vector x(0) arbitrarily, and compute x(1) from x(1) = Tjx(0) + cj. We then use x(1) to update x(0), and the process is repeated.

Before executing the program, we store the matrix A into the matrix vaiable [A]. The vector b is stored into the n by 1 "column matrix" [B], and the initial iterate x(0) is stored into the n by 1 column matrix [C]. The program will compute the matrix D and store it into [D] (see instructions#3,4,5,6,7,8,9). The matrix D - A is the same as L + U. This matrix L + U is stored into [E] (see instruction#10). The matrix Tj is stored into [F], and the vector cj is stored into [G] (see instructions#11,12).

Since we need to compute the matrix Tj and the vector cj only once, we use the entry [A](1,1) as a signal for the computation of the splitting. Initially, [A](1,1) is not zero, as required by the Jacobi method. After the splitting is computed, we save [A](1,1) into the variable S (see instruction#13) and set [A](1,1) to zero (see instruction#14). Thus, if [A](1,1) is zero, the program will skip the computation of the splitting (see instructions#1,2,15), since this computation has already been done.

The program will next compute x(1) and store x(1) into [H] (see instruction#16). The difference vector x(1) - x(0) is stored into [I] (see instruction#17). x(1) is used to update x(0) (see instruction#18). An output matrix [J] is constructed in instruction#19. The first column of [J] is the vector x(1). The second column of [J] is the previously computed vector x(1) - x(0). Instruction#19 also displays the output matrix [J].

After executing the program, keep pressing [ENTER]. The iterates will appear on the home screen. In general, convergence is linear.

(Note: The above matrix description of Jacobi iteration is, in general, suitable for solving small size textbook problems. For large sparse systems, it is generally more efficient computationally to use the "equation by equation" algorithm for Jacobi iteration, as described in most textbooks on numerical analysis.)

Instruction#1_____If [A](1,1) = 0..........(See note at.end of program.)
Instruction#2_____Go to 1
Instruction#3_____dim([A]) -> L1..........(Note: We use the symbol "->" to
Instruction#4_____L1 -> dim([D])...........represent the operation "Store".)
Instruction#5_____Fill(0,[D])
Instruction#6_____L1(1) -> N
Instruction#7_____For(I,1,N)
Instruction#8_____[A](I,I) -> [D](I,I)
Instruction#9_____End
Instruction#10____[D] - [A] -> [E]
Instruction#11____[D]-1[E] -> [F]
Instruction#12____[D]-1[B] -> [G]
Instruction#13____[A](1,1) -> S
Instruction#14____0 -> [A](1,1)
Instruction#15____Lbl 1
Instruction#16____[F][C] + [G] -> [H]
Instruction#17____[H] - [C] -> [I]
Instruction#18____[H] -> [C]
Instruction#19____augment([H],[I]) -> [J]

Note: Use the MATRIX menu to enter the matrix names [A] , [B] , ... onto a program line. Do not type the characters "[" , "A" , "]" , "B" , ... separately.

-------------------------------------------------------------------------

Program #12 GSEIDEL (Gauss-Seidel Iteration)

Gauss-Seidel iteration (Seidel iteration for short) is a method for solving a linear system of equations Ax = b , where A is an n by n matrix. In matrix form, Seidel iteration may be described as follows:

We "split" the matrix A into D - L - U, where D is a diagonal matrix with Dii = Aii, L is a strictly lower triangular matrix with Lij = -Aij when i > j, and U is a strictly upper triangular matrix with Uij = -Aij when j > i. We require that Aii is not zero for each i. Then, (D-L)-1 exists, and AX = b may be rewritten as x = (D-L)-1U x + (D-L)-1b. If we set Tg = (D-L)-1U and cg = (D-L)-1 b , then the original system is equivalent to the system x = Tgx + cg . If the spectral radius of Tg is less than 1, then we can solve by iteration. Thus, we choose the vector x(0) arbitrarily, and compute x(1) from x(1) = Tgx(0) + cg. We then use x(1) to update x(0), and the process is repeated.

Before executing the program, we store the matrix A into the matrix vaiable [A]. The vector b is stored into the n by 1 "column matrix" [B], and the initial iterate x(0) is stored into the n by 1 column matrix [C]. The program will compute the matrix D-L and store it into [D] (see instructions#3,4,5,6,7,8,9,10,11). The matrix D -L - A is the same as U and is stored into [E] (see instruction#12). The matrix Tg is stored into [F], and the vector cg is stored into [G] (see instructions#13,14).

Since we need to compute the matrix Tg and the vector cg only once, we use the entry [A](1,1) as a signal for the computation of the splitting. Initially, [A](1,1) is not zero, as required by the Seidel method. After the splitting is computed, we save [A](1,1) into the variable S (see instruction#15) and set [A](1,1) to zero (see instruction#16). Thus, if [A](1,1) is zero, the program will skip the computation of the splitting (see instructions#1,2,17), since this computation has already been done.

The program will next compute x(1) and store x(1) into [H] (see instruction#18). The difference vector x(1) - x(0) is stored into [I] (see instruction#19). x(1) is used to update x(0) (see instruction#20). An output matrix [J] is constructed in instruction#21. The first column of [J] is the vector x(1). The second column of [J] is the previously computed vector x(1) - x(0). Instruction#21 also displays the output matrix [J].

After executing the program, keep pressing [ENTER]. The iterates will appear on the home screen. In general, convergence is linear.

(Note: The above matrix description of Seidel iteration is, in general, suitable for solving small size textbook problems. For large sparse systems, it is generally more efficient computationally to use the "equation by equation" algorithm for Seidel iteration, as described in most textbooks on numerical analysis.)

Instruction#1_____If [A](1,1) = 0.............(See note at end of program.)
Instruction#2_____Go to 1
Instruction#3_____dim([A]) -> L1..........(Note: We use the symbol "->" to
Instruction#4_____L1 -> dim([D])...........represent the operation "Store".)
Instruction#5_____Fill(0,[D])
Instruction#6_____L1(1) -> N
Instruction#7_____For(I,1,N)
Instruction#8_____For(J,1,I)
Instruction#9_____[A](I,J) -> [D](I,J)
Instruction#10____End
Instruction#11____End
Instruction#12____[D] - [A] -> [E]
Instruction#13____[D]-1[E] -> [F]
Instruction#14____[D]-1[B] -> [G]
Instruction#15____[A](1,1) -> S
Instruction#16____0 -> [A](1,1)
Instruction#17____Lbl 1
Instruction#18____[F][C] + [G] -> [H]
Instruction#19____[H] - [C] -> [I]
Instruction#20____[H] -> [C]
Instruction#21____augment([H],[I]) -> [J]

Note: Use the MATRIX menu to enter the matrix names [A] , [B] , ... onto a program line. Do not type the characters "[" , "A" , "]" , "B" , ... separately.

-------------------------------------------------------------------------

Program #13 SOR (Successive Over-Relaxation)

Successive over-relaxation iteration (SOR iteration for short) is a method for solving a linear system of equations Ax = b , where A is an n by n matrix. In matrix form, SOR iteration may be described as follows:

We "split" the matrix A into D - L - U, where D is a diagonal matrix with Dii = Aii, L is a strictly lower triangular matrix with Lij = -Aij when i > j, and U is a strictly upper triangular matrix with Uij = -Aij when j > i. We require that Aii is not zero for each i. Then, (D-wL)-1 exists, where the parameter w is an arbitrary number. The system AX = b may then be rewritten as x = (D-wL)-1[(1 - w) D + w U] x + w(D-wL)-1b. If we set Tw = (D-wL)-1[(1 - w) D + w U] and cw = w(D-wL)-1 b , then the original system is equivalent to the system x = Twx + cw . If the spectral radius of Tw is less than 1, then we can solve by iteration. Thus, we choose the vector x(0) arbitrarily, and compute x(1) from x(1) = Twx(0) + cw. We then use x(1) to update x(0), and the process is repeated.

Before executing the program, we store the matrix A into the matrix vaiable [A]. The vector b is stored into the n by 1 "column matrix" [B], and the initial iterate x(0) is stored into the n by 1 column matrix [C]. The parameter w is stored into the variable W. The program will first compute the matrices D and D-L and store them into the matrix variables [D] and [E] respectively (see instructions#3,4,5,6,7,8,9,10,11,12,13. Instructions#3,4,5,6 initially fill [D] and [E] with zeros.). The matrix D -L - A is the same as U and is stored into [F] (see instruction#14). The matrix D - (D - L) is the same as L and is stored back into [E] (see instruction#15). At this point, [D], [E] and [F] contain the matrices D, L and U respectively. The matrix Tw is stored into [G], and the vector cw is stored into [H] (see instructions#16,17,18).

Since we need to compute the matrix Tw and the vector cw only once, we use the entry [A](1,1) as a signal for the computation of the splitting. Initially, [A](1,1) is not zero, as required by the SOR method. After the splitting is computed, we save [A](1,1) into the variable S (see instruction#19) and set [A](1,1) to zero (see instruction#20). Thus, if [A](1,1) is zero, the program will skip the computation of the splitting (see instructions#1,2,21), since this computation has already been done.

The program will next compute x(1) and store x(1) into [I] (see instruction#22). The difference vector x(1) - x(0) is stored into [J] (see instruction#23). x(1) is used to update x(0) (see instruction#24). The matrices [I] and [J] are used to construct an output matrix, and the result is put back into [J] (see instruction#25). Thus, the first column of the output matrix [J] is the vector x(1), and the second column of [J] is the previously computed vector x(1) - x(0). Instruction#25 also displays the output matrix [J].

After executing the program, keep pressing [ENTER]. The iterates will appear on the home screen. In general, convergence is linear.

(Note: The above matrix description of SOR iteration is, in general, suitable for solving small size textbook problems. For large sparse systems, it is generally more efficient computationally to use the "equation by equation" algorithm for SOR iteration, as described in most textbooks on numerical analysis.)

Instruction#1_____If [A](1,1) = 0...........(See note at end of program.)
Instruction#2_____Go to 1
Instruction#3_____dim([A]) -> L1..........(Note: We use the symbol "->" to
Instruction#4_____L1 -> dim([D])...........represent the operation "Store".)
Instruction#5_____Fill(0,[D])
Instruction#6_____[D] -> [E]
Instruction#7_____L1(1) -> N
Instruction#8_____For(I,1,N)
Instruction#9_____[A](I,I) -> [D](I,I)
Instruction#10____For(J,1,I)
Instruction#11____[A](I,J) -> [E](I,J)
Instruction#12____End
Instruction#13____End
Instruction#14____[E] - [A] -> [F]
Instruction#15____[D] - [E] -> [E]
Instruction#16____( [D] - W [E] )-1 -> [J]
Instruction#17____( [J] )( ( 1-W ) [D] + W [F] ) -> [G]
Instruction#18____( W [J] )( [B] ) -> [H]
Instruction#19____[A](1,1) -> S
Instruction#20____0 -> [A](1,1)
Instruction#21____Lbl 1
Instruction#22____[G][C] + [H] -> [I]
Instruction#23____[I] - [C] -> [J]
Instruction#24____[I] -> [C]
Instruction#25____augment([I],[J]) -> [J]

Note: Use the MATRIX menu to enter the matrix names [A] , [B] , ... onto a program line. Do not type the characters "[" , "A" , "]" , "B" , ... separately.

-------------------------------------------------------------------------

Program #14 SCLPWR (Scaled Power Method)

Let A be an n by n real matrix. We assume there is a unique dominant eigenvalue k1, so that abs(k1) > abs(k2) > ... > abs(kn) . In the scaled power method, we start with a vector x(0) , which we assume to be an approximation to the eigenvector v(1) corresponding to the eigenvalue k1. (In practice, x(0) is usually chosen arbitrarily. However, if an approximation to the eigenvector v(1) is available, this approximation should be used for x(0).) The vector x(0) is scaled as follows:

Let xi(0) , i=1, ... , n be the components of x(0). Then, the set of numbers abs(xi(0)) has a maximum value. Let p be the smallest index such abs(xp(0)) is this maximum value. Then, to scale the vector x(0), we multiply x(0) by (xp(0))-1 . Then, the scaled vector has components between -1 and +1, and at least one component is equal to +1. Also, the pth component of the scaled vector is +1.

Now assume that x(0) is already scaled, where xp(0) = +1. We form the product Ax(0) and define this vector to be x(1). If x(0) approximates the eigenvector v(1), then x(1) approximates the vector k1v(1). Since xp(0) = 1, then vp(1) is approximately 1. Thus, the pth component of k1v(1) is approximately k1. Thus, xp(1) is approximately the dominant eigenvalue k1. At this point, we replace the vector x(0) with x(1), and the process is repeated.

Before executing the program, we store the matrix A into the matrix variable [A]. The vector x(0) is stored in [B]. The dimension of A is determined in instructions#1,2. The absolute values of the components of x(0) are stored in [C] (see instruction#3). We determine the index p in instructions#4,5,6,7,8,9,10,11,12. This index p is stored in the variable P. The vector x(0) is scaled (see instructions#13,14), and the scaled vector is stored back into [B]. The scaled vector is also stored in the list L2 (see instruction#15) in order to prepare for the output of the program.

The vector x(1) = Ax(0) is computed and stored back into [B] and also into the list L3 (see instructions#16,17). A list of zeroes of length n is created (see instructions#18,19), and the pth component of this list is then replaced by xp(1), the approximation to the dominant eigenvalue k1 (see instruction#20). An output matrix [C] is set up in instruction#21. In this matrix, the first column is the scaled vector x(0), the second column is the vector x(1), and the third column contains xp(1), the approximation to the eigenvalue k1. The output matrix [C] is displayed in instruction#22.

After executing the program, keep pressing [ENTER]. The iterates for x(0), x(1), and approximations to k1 appear on the home screen. Convergence to k1 is, in general, linear.

Instruction#1_____dim([A]) -> L1..........(See note at end of program.)
Instruction#2_____L1(1) -> N..................(Note: We use the symbol "->" to
Instruction#3_____abs([B]) -> [C].............represent the operation "Store".)
Instruction#4_____[C](1,1) -> X
Instruction#5_____1 -> P
Instruction#6_____For ( I , 2 , N )
Instruction#7_____If ( [C](I,1) - X ) > 0
Instruction#8_____Then
Instruction#9_____I -> P
Instruction#10____[C](P,1) -> X
Instruction#11____End
Instruction#12____End
Instruction#13____[B](P,1) -> X
Instruction#14____X-1[B] -> [B]
Instruction#15____Matr>list([B],L2)
Instruction#16____[A][B] -> [B]
Instruction#17____Matr>list([B],L3)
Instruction#18____N -> dim(L4)
Instruction#19____Fill(0,L4)
Instruction#20____[B](P,1) -> L4(P)
Instruction#21____List>matr(L2,L3,L4,[C])
Instruction#22____[C]

Note: Use the MATRIX menu to enter the matrix names [A] , [B] , ... onto a program line. Do not type the characters "[" , "A" , "]" , "B" , ... separately.

-------------------------------------------------------------------------

Program #15 WDEFLATE (Wielandt Deflation)

We assume A is an n by n matrix with distinct eigenvalues ki , i = 1, ... , n , and corresponding eigenvectors v(i) , i = 1, ... , n. We also assume abs(k1) > abs(k2) > ... , abs(kn) . Let x be a vector with the property that xTv(1) = 1. Let B = A - k1v(1)xT. Then it is not hard to show that B has eigenvalues 0, k2, ... , kn . Also, v(1) is an eigenvector of B corresponding to the eigenvalue 0 . Also, if w(2), ... , w(n) are eigenvectors of B corresponding to k2, ... , kn, then it can be shown that v(i) = ( ki - k1 ) w(i) + k1 (xT w(i) ) v(1) , i = 2 , ... , n .

We make the simplifying assumption for Wielandt deflation that the nth component of v(1) is not zero. That is, we assume vn(1) is not zero. We also assume k1 is not zero. Then we use the last row of the matrix A to help define the vector x. We set x = ( k1 vn(1) )-1 * [ an1, an2, ... , ann ]T . Thus, x is a multiple of the last row of A. With this choice of x, it is not hard to show that xT v(1) = 1. Also, it is not hard to show that the last row of B is zero. As a consequence, the last component of each eigenvector w(2), ... , w(n) of B must be zero.

As a resut of the above, the eigenvalues k2, ... , kn of B are the eigenvalues of B' , the upper left (n-1) by (n-1) block of B. If w' is an eigenvector of B', then w = [ w' , 0 ]T is an eigenvector of B.

(Thus, if we know k1 and v(1) (from a method such as the power method) for the matrix A, we can reduce the problem of finding k2 and v(2) for A as follows: We consider the smaller size matrix B', we find k2 for B' (from some method such as the power method), and then we find the corresponding eigenvector ( w(2) )'. With this eigenvector, we easily construct w(2) and v(2). The process of finding B' is called Wielandt deflation. The process of deflation could be continued to find B'', etc... .)

Before executing our program, we put the matrix A into the matrix variable [A]. We also put the vector v(1) into the column matrix [C]. Instructions#1,2 store the dimension n of [A] into the variable N. Instructions#3,4,5,6 put the last row of [A] into the column matrix [E]. Instructions#7,8 use the equation [A v(1)]n = [k1 v(1)]n to compute the value of k1 and store this value in the variable K. Instruction#9 forms the vector x. Instructions#10,11 form the matrix B. (Instructions#7,8,9,10,11 could be modified to save some computation.) Instruction#12 reduces the size of B, so that the last row and column of B are deleted. Thus, B' is left in the matrix variable [B]. Instruction#13 displays B'.

Instruction#1_____dim([A]) -> L1..........(See note at end of program.)
Instruction#2_____L1(1) -> N.........................(Note: We use the symbol "->" to
Instruction#3_____dim([C]) -> dim([E])..........represent the operation "Store".)
Instruction#4_____For (I,1,N)
Instruction#5_____[A](N,I) -> [E](I,1)
Instruction#6_____End
Instruction#7_____[E]T[C] -> [D]
Instruction#8_____[D](1,1)/[C](N,1) -> K
Instruction#9_____( [D](1,1) )-1 [E] -> [E]
Instruction#10____[C][E]T -> [F]
Instruction#11____[A] - K [F] -> [B]
Instruction#12____{N-1,N-1} -> dim( [B] )
Instruction#13____[B]

Note: Use the MATRIX menu to enter the matrix names [A] , [B] , ... onto a program line. Do not type the characters "[" , "A" , "]" , "B" , ... separately.

-------------------------------------------------------------------------

Program #16 INTERPN (N-Point Interpolation)

We assume the n points are ( x0 , f0 ) , ... , (xn-1 , fn-1 ) . Then, there is a unique interpolating polynomial P(x) = a0 + a1x + ... + an-1xn-1 which passes through the n data points. This program computes the values of a0, ... , an-1 .

Before executing the program, put the list { x0 , x1 , ... , xn-1 } into the list variable L1, and { f0 , f1 , ... , fn-1 } into the list variable L2. Thus, L1 and L2 should each contain only n numbers. Instructions #1,2,3 set up matrices [A] and [B] to be used by the program. Instructions #4,5,6,7,8,9,10 set up an augmented Vandermonde matrix in [A]. Thus, Aij = xi-1^(j-1) , and the last column of [A] contains f0, ... , fn-1 . Instruction #11 reduces [A] to its row reduced echelon form, and thus solves the system P(xi) = fi , i = 0, 1, ... , n-1 for a0, ... , an-1 . (Note: The instruction rref may be found in the MATRIX, MATH menu.) The solution is left in the last column of [A]. (Since P(x) exists and is unique, [A] must be invertible.)

Instructions #12,13,14 move the last column of [A] into the column matrix [B]. Thus, [B] contains the column vector (a0 , ... , an-1 )T . Instruction#15 displays [B] on the home screen. (Note: If one wants to reverse the order of the numbers in [B], then instruction#13 should be replaced with [A](I,N+1) -> [B](N-I+1,1).)

When the number of data points is n = 2,3,4 or 5, then the TI-83 regression functions make it possible to write programs for interpolation which are very short. See the programs INTERP2, INTERP3, INTERP4, and INTERP5.

Instruction#1_____dim(L1) -> N..................(Note: We use the symbol "->" to
Instruction#2_____{N,N+1} -> dim([A]).......represent the operation "Store".)
Instruction#3_____{N,1} -> dim([B])............(See note below.)
Instruction#4_____For (I,1,N)
Instruction#5_____L2(I) -> [A](I,N+1)
Instruction#6_____1 -> [A](I,1)
Instruction#7_____For (J,2,N)
Instruction#8_____L1(I)^(J-1) -> [A](I,J)
Instruction#9_____End
Instruction#10____End
Instruction#11____rref([A]) -> [A]
Instruction#12____For (I,1,N)
Instruction#13____[A](I,N+1) -> [B](I,1)
Instruction#14____End
Instruction#15____[B]

Note: Use the MATRIX menu to enter the matrix names [A] , [B] , ... onto a program line. Do not type the characters "[" , "A" , "]" , "B" , ... separately.

-------------------------------------------------------------------------

Program #17 INTERP2 (2-Point Interpolation)

For INTERP2, we assume the values for x0, x1 are stored in the list variable L1, and the values for f0, f1 are stored in the list variable L2. We also assume that a value for x is stored in the variable X. Thus, L1 and L2 should each contain only two numbers. (Before setting up the lists L1,L2, clear these lists in order to be sure the numbers you enter are treated as real numbers. Also, clear any expression in the Y1 function variable.) Instruction#1 computes the linear interpolating polynomial y = a + bx. (In instruction#1, the function LinReg(a+bx) may be found in the STAT, CALC menu, item#8. We could equally well use item#4, LinReg(ax+b).) The polynomial a + bx is stored in the function variable Y1 when instruction#1 is executed. Instruction#2 evaluates Y1 at the given value of X and displays the result on the home screen. The values for a,b may be read from Y1 in the "Y=" menu. (Note: Set mode to display 10 decimal places before executing the program.) Or, you can display the contents of "a" and "b" on the home screen by selecting "a" and "b" in the VARS menu.

Instruction#1_____LinReg(a+bx)L1,L2,Y1
Instruction#2_____Y1

-------------------------------------------------------------------------

Program #18 INTERP3 (3-Point Interpolation)

The description for INTERP3 is very much the same as the description for INTERP2. Only minor changes occur in the description for INTERP3.

For INTERP3, we assume the values for x0, x1, x2 are stored in the list variable L1, and the values for f0, f1, f2 are stored in the list variable L2. We also assume that a value for x is stored in the variable X. Thus, L1 and L2 should each contain only three numbers. (Before setting up the lists L1,L2, clear these lists in order to be sure the numbers you enter are treated as real numbers. Also, clear any expression in the Y1 function variable.) Instruction#1 computes the quadratic interpolating polynomial y = ax2 + bx + c . (In instruction#1, the function QuadReg may be found in the STAT, CALC menu, item#5.) The polynomial ax2 + bx + c is stored in the function variable Y1 when instruction#1 is executed. Instruction#2 evaluates Y1 at the given value of X and displays the result on the home screen. The values for a,b,c may be read from Y1 in the "Y=" menu. (Note: Set mode to display 10 decimal places before executing the program.) Or, you can display the contents of a,b,c on the home screen by selecting a,b,c in the VARS menu.

Instruction#2_____Y1

-------------------------------------------------------------------------

Program #19 INTERP4 (4-Point Interpolation)

The description for INTERP4 is very much the same as the description for INTERP3. Only the obvious minor changes need to be made.

The program for INTERP4 is as follows:

Instruction#1_____CubicReg L1,L2,Y1
Instruction#2_____Y1

-------------------------------------------------------------------------

Program #20 INTERP5 (5-Point Interpolation)

The description for INTERP5 is very much the same as the description for INTERP3. Only the obvious minor changes need to be made.

The program for INTERP5 is as follows:

Instruction#1_____QuartReg L1,L2,Y1
Instruction#2_____Y1

-------------------------------------------------------------------------

Program #21 NEVILLE (Neville's Method)

In Neville's method, we assume that we have n data points (x0,f0), ... , (xn-1,fn-1). If a value for x is given, then Neville's method computes the value of the interpolating polynomial P(x) which passes through the n points.

Nevile's method is based on Aitken's Lemma. Aitken's Lemma may be described as follows: Let V be the set of indices {0,1,...,n-1} of the x-coordinates of the data points. Let {i,j,...,k} be a subset of V. Define Pij...k(x) ( =PV(x) ) to be the interpolating polynomial which uses the data points at xi, xj, ... , xk . Let S and T be subsets of V with the following properties:

1. S and T have the same number of indices.
2. Every index in S is also in T except for one index. Call this index "i" (in S).

Then it follows that every index in T is also in S except for one index. Call this index "j" (in T). It also follows that "i" is not equal to "j".

The main statement for Aitken's Lemma is the following:

Let U be the union of S and T. Then

PU(x) = [ ( x-xj )PS(x) - ( x-xi )PT(x) ] / [ (x-xj ) - ( x-xi ) ] .

We call this equation "Aitken's formula". (The denominator in this equation can be simplified to xi - xj (a constant) if so desired.)

Neville's method makes use of Aitken's Lemma a total of n(n-1)/2 times.The computation is organized as follows: Let [A] be an n by n matrix. All entries above the diagonal are set to zero. In the first column, we place the numbers f0, f1, ... , fn-1 . These numbers are also the values of P0(x), P1(x), ... , Pn-1(x) (for the given value of x. Remember, Pk(x) is a zeroth degree polynomial using our notation.) We now use P0, P1, ... , Pn-1 pairwise in Aitken's Lemma to compute P01, P12, ... , Pn-2,n-1 (for the given value of x. For example, P01 = [ ( x-x1 )*P0 - ( x-x0 )*P1 ] / [ ( x-x1 ) - ( x-x0 ) ] . ) These numbers are placed into the 2nd column of [A] , on and below the diagonal. These numbers in the second column column are now used pairwise in Aitken's Lemma to compute P012, P123, ... , Pn-3,n-2,n-1 . (For example, P012 = [ ( x-x2 )*P01 - ( x-x0 )*P12 ] / [ ( x-x2 ) - ( x-x0 ) ] . ) These new numbers are placed into the 3rd column of [A], on and below the diagonal. This process continues until we reach the lower right corner of [A]. This entry will contain P0,1,2,...,n-1 , the value of the interpolating polynomial using the complete set of data.

In the following program, we assume x0, ... , xn-1 are in the list variable L1, and f0, ... , fn-1 are in the list variable L2. We also assume the value of x is in the variable X. Instruction#1 finds the number of data points n and puts the result into N. Instruction#2 sets up an (n) by (n+2) matrix [A] which will be used as an output matrix. Instruction#3 sets the initial values of the entries of [A] to zero. Instructions#4,5,6,7,8 form a loop which puts the numbers {xi} into the (n+2)nd column of [A], puts the numbers {x-xi} (for use in Aitken's formula) into the (n+1)st column of [A], and puts {fi} into the 1st column of [A].

Instructions#9,10,11,12,13 form a double loop which fills columns 2,3,...,n of [A] according to Neville's method as described above. Instuction#11 is the main instruction. This long instruction uses Aitken's formula as follows: Two consecutive numbers are taken from column (j-1) of [A]. These numbers are in rows (i-1) and (i) of [A]. Also, these numbers are PS(x) and PT(x) in Aitken's formula. These numbers are used with the two appropriate values from column (n+1) of [A] (which contains the values of {x-xk} ).The result from Aitken's formula is stored into the (i,j) entry of [A].

Instruction#14 displays the entire matrix [A] on the home screen. The (n,n) entry of [A] contains the value of the interpolating polynomial which uses the entire set of data points.

Before executing the program, put {x0, ... , xn-1} into L1, {f0, ... , fn-1} into L2, and x into X. Then execute the program and read the results on the home screen. You will probably need to use the arrow keys to scroll through the matrix. Or, you can go to the matrix editor and read the entries there. (In the matrix editor, you can read extra decimal places at the bottom of the screen.)

Instruction#1_____dim(L1) -> N..................(Note: We use the symbol "->" to
Instruction#2_____{N,N+2} -> dim([A]).......represent the operation "Store".)
Instruction#3_____Fill(0,[A]).........................(See note below.)
Instruction#4_____For (I,1,N)
Instruction#5_____L1(I) -> [A](I,N+2)
Instruction#6_____X-L1(I) -> [A](I,N+1)
Instruction#7_____L2(I) -> [A](I,1)
Instruction#8_____End
Instruction#9_____For (J,2,N)
Instruction#10____For (I,J,N)
Instruction#11____( [A](I-1,J-1)*[A](I,N+1) - [A](I,J-1)*[A](I-J+1,N+1) ) /
________________( [A](I,N+1) - [A](I-J+1,N+1) ) -> [A](I,J)
Instruction#12____End
Instruction#13____End
Instruction#14____[A]

Note: Use the MATRIX menu to enter the matrix name [A] onto a program line. Do not type the characters "[" , "A" , "]" separately.

-------------------------------------------------------------------------

Program #22 NEWTDIV (Newton's Divided Difference Method)

In Newton's divided difference method, we assume that we have n data points (x0,f0), ... , (xn-1,fn-1). Then there is a unique interpolating polynomial P(x) of degree (n-1) which passes through the n data points. This polynomial may be expressed in the form a0 + a1*(x-x0) + a2*(x-x0)*(x-x1) + ... + an-1*(x-x0)*(x-x1)* ... *(x-xn-2) . The constants a0, ... , an-1 are given by ak = f [ x0, x1, ... , xk-1 ] , where f [ x0, x1, ... , xk-1 ] is the Newton divided difference of order (k-1) using the points x0, x1, ... , xk-1 . This program computes the values of a0, ... , an-1.

The Newton divided differences are defined recursively. The order-0 divided differences are defined by f[xj] = fj, j=0, 1, ... , n-1. Thus, there are n divided differences of order-0. These n numbers will be placed into the first column of an n by n matrix [A]. (The entries above the diagonal of [A] are set to zero.) The order-k divided differences are defined by

f[ xi, xi+1, ... , xi+k ] = ( f[ xi+1, xi+2 , ... , xi+k ] - f[ xi, xi+1, ... , xi+k-1 ] ) / ( xi+k - xi ) for i = 0, 1, ... , n - k - 1 .

Thus, there are n-k such numbers. These will be placed into column# (k+1) of the matrix [A] , on and below the diagonal.

Thus, column#2 will contain the numbers f[x0,x1], f[x1,x2], ... , f[xn-2, xn-1] , where f[x0,x1] = (f[x1] - f[x0]) / (x1-x0), f[x1,x2] = (f[x2] - f[x1]) / (x2-x1) , etc... .

Column#3 will contain the numbers f[x0,x1,x2], f[x1,x2,x3], ... , f[xn-3,xn-2, xn-1] , where f[x0,x1,x2] = (f[x1,x2] - f[x0,x1]) / (x2-x0), f[x1,x2,x3] = (f[x2,x3] - f[x1,x2]) / (x3-x1) , etc... .

Column#n of [A] will contain the single number f [ x0, x1, ... , xn-1 ] in the diagonal entry [A]n,n .

Before executing the program, put the list {x0, ... , xn-1} into the list variable L1, {f0, ... , fn-1} into the list variable L2, and a value for x into the variable X. Instruction#1 finds the number of data points n and puts the result into N. Instruction#2 sets up an (n) by (n+2) matrix [A] which will be used as an output matrix. Instruction#3 sets the initial values of the entries of [A] to zero. Instructions#4,5,6,7 form a loop which puts the numbers {xi} into the (n+1)st column of [A], and puts {fi} into the 1st column of [A].

Instructions#8,9,10,11,12 form a double loop which fills columns 2, 3, ... , n of [A] with divided differences of order 1, 2, ... , (n-1) respectively (on and below the diagonal). Instruction#10 is the main instruction. This long instruction uses two consecutive divided differences in column (j-1) of [A] to form a next-order divided difference to be placed into column (j) of [A]. After the double loop is completed, the numbers a0, a1, ... , an-1 will appear on the main diagonal of [A].

Instructions#13,14,15,16,17,18 are optional. They compute the value of P(x) at the given value of x using a Horner type nest. Thus, P(x) is written in the form

P(x) = ( ... (( an-1*(x-xn-1) + an-2 )*(x-xn-2) + an-3 )*(x-xn-3) + ... )*(x-x0) + a0 .

The value of the Horner sum is P(x), and this number is accumulated in the variable S (instruction#16) and also placed into the first entry of column (n+2) of [A] (instruction#18).

Instruction#19 displays the entire matrix [A] on the home screen. You will probably need to use the arrow keys to scroll through the matrix. Or, you can go to the matrix editor and read the matrix entries there. (In the matrix editor, you can read extra decimal places at the bottom of the screen.)

Instruction#1_____dim(L1) -> N..................(Note: We use the symbol "->" to
Instruction#2_____{N,N+2} -> dim( [A] ).......represent the operation "Store".)
Instruction#3_____Fill( 0, [A] )......................(See note below.)
Instruction#4_____For(I,1,N)
Instruction#5_____L1(I) -> [A](I,N+1)
Instruction#6_____L2(I) -> [A](I,1)
Instruction#7_____End
Instruction#8_____For(J,2,N)
Instruction#9_____For(I,J,N)
Instruction#10____( [A](I,J-1) - [A](I-1,J-1) ) / ( L1(I) - L1(I-J+1) ) -> [A](I,J)
Instruction#11____End
Instruction#12____End
Instruction#13____[A](N,N) -> S
Instruction#14____For(I,2,N)
Instruction#15____N-I+1 -> K
Instruction#16____S * ( X - L1(K) ) + [A](K,K) -> S
Instruction#17____End
Instruction#18____S -> [A](1,N+2)
Instruction#19____[A]

Note: Use the MATRIX menu to enter the matrix name [A] onto a program line. Do not type the characters "[" , "A" , "]" separately.

-------------------------------------------------------------------------

Program #23 LAGRANGE (Lagrange Interpolation Polynomial)

In Lagrange's method, we assume that we have n data points (x0,f0), ... , (xn-1,fn-1). Then there is a unique interpolating polynomial P(x) of degree (n-1) which passes through the n data points. This polynomial may be expressed in the form P(x) = f0L0(x) + f1L1(x) + ... + fn-1Ln-1(x), where

Li(x) = [ (x-x0) (x-x1)...(x-xi-1) (x-xi+1)...(x-xn-1)] / [ (xi-x0) (xi-x1)...(xi-xi-1) (xi-xi+1)...(xi-xn-1) ] .

Notice that in the numerator of Li(x), the factor (x-xi) does not appear. Also, the denominator is obtained from the numerator by replacing x with xi. ( Li(x) is called the ith Lagrange interpolating coefficient polynomial. This coefficient polynomial multiplies fi when forming P(x). Li(x) has the property that Li(xj) = 0 if "i" is not equal to "j", and Li(xi) = 1 .)

For a given value of x (not equal to any of the xi values), this program will compute a list of the values of the numerators of the functions Li(x), a list of the values of the denominators of the functions Li(x), a list of the values of Li(x), a list of the values of fiLi(x), and a value for P(x).

Before executing the program, put the list {x0, ... , xn-1} into the list variable L1, {f0, ... , fn-1} into the list variable L2, and a value for x (not equal to any of the xi values) into the variable X. Instruction#1 finds the number of data points n and puts the result into N. Instruction#2 creates a list L4 with n entries, which will be used by the program to store the denominators of Li(x), i = 0, ... , n-1 . Instruction#3 forms the product (x-x0) (x-x1) ... (x-xn-1) and puts the result into Q. Q will be used to help form the numerators of Li(x) .

Instructions#4,5,6,7,8,9 are optional. They test to see if Q = 0 (instruction#4), which would mean that x = xi for some i. Since this would cause a division by zero in instruction#10, a message is displayed to change the value of x (instructions#6,7). Instructions #4,5,6,7,8,9 could be deleted. Then, if x = xi for some i, the program would stop at instruction#10 because of a division by zero.

In instruction#10, the term (x-L1) in the denominator creates a list { x-x0 , x-x1 , ... , x-xn-1 } . Division of Q by this list creates another list { Q/(x-x0) , Q/(x-x1) , ... , Q/(x-xn-1) } . Since Q is the product (x-x0) (x-x1) ... (x-xn-1) , this last list contains the numerators of the Lagrange coefficients L0(x), ... , Ln-1(x) . This list of numerators is stored into the list variable L3.

Instructions#11,12,13,14,15,16,17 form a double loop which computes the denominators of the Lagrange coefficients Li(x) and stores them into the list variable L4. In the inner loop (instructions #13,14,15,16), J is tested to see if J = I (instruction#14). If J = I, the factor (xi-xi) (which equals zero) is left out of the denominator of Li(x), since instruction#15 is skipped when J = I .

Instruction#18 computes the values of Li(x) and puts them into the list variable L5. Instruction#19 computes the values of fiLi(x) and puts them into the list variable L6. Instruction#20 sums the values of fiLi(x), which results in the value of P(x), and stores this result in P. Instruction#21 sets up an output matrix in the matrix variable [A]. The first column of [A] contains the value of P(x) (followed by zeroes). Column#2 contains {xi}. Column#3 contains {fi}. Column#4 contains the numerators of Li(x). Column#5 contains the denominators of Li(x). Column#6 contains the values of Li(x). The last column (column#7) contains the values of fiLi(x).

Instruction#22 displays the entire matrix [A] on the home screen. You will probably need to use the arrow keys to scroll through the matrix. Or, you can go to the matrix editor and read the matrix entries there. (In the matrix editor, you can read extra decimal places at the bottom of the screen.)

Instruction#1_____dim(L1) -> N.............(Note: We use the symbol "->" to
Instruction#2_____N -> dim(L4).............represent the operation "Store".)
Instruction#3_____prod(X-L1) -> Q
Instruction#4_____If Q = 0
Instruction#5_____Then
Instruction#6_____Disp "X=X(I) FOR SOME"
Instruction#7_____Disp "I. CHANGE X"
Instruction#8_____Stop
Instruction#9_____End
Instruction#10____Q/(X-L1) -> L3
Instruction#11____For (I,1,N)
Instruction#12____1 -> L4(I)
Instruction#13____For (J,1,N)
Instruction#14____If J != I ________________(See note below)
Instruction#15____L4(I)*(L1(I) - L1(J) ) -> L4(I)
Instruction#16____End
Instruction#17____End
Instruction#18____L3 / L4 -> L5
Instruction#19____L2L5 -> L6
Instruction#20____sum(L6) -> P
Instruction#21____List>matr( {P} , L1, L2, L3, L4, L5, L6, [A] )....(See note below.)
Instruction#22____[A]

Note that in instruction#14, we are using the symbol != to represent "not equal". When entering this program into the TI83, one should use the more standard symbol from the TEST menu.
For instructions#21,22, use the MATRIX menu to enter the matrix name [A] onto a program line. Do not type the characters "[" , "A" , "]" separately.