Math 438/538 Numerical Analysis Spring 2011: Brian Hassard

Course Materials:

Classes are Tuesday-Thursday, 12:30-1:50 in Math 150. The first class is Tuesday, January 18.
Labs (recitations) start Wednesday, January 26 and will be Wednesdays, 8:00-8:50 in Math 250.

To review the material from MTH 537 (Intro Numerical Analysis 1) Fall 2010, please see the website for last term (MTH 437-537) is http://www.math.buffalo.edu/~hassard/437-537/

Program eulers_method_with_exact_and_err.m for y' = -2 y, y(0) = 1 on [a,b]=[0,1]

Output from running program (octave) for stepsizes h=0.1, 0.01, 0.001, 0.0001, 0.00001

The code euler2.m applies Euler's method to a system of 2 first order ODEs.
The code is modular, consisting of function euler2, function eulerstep, and function ydot.
Within euler2, n is the number of steps of size h,
the values t(1) through t(n+1) correspond to
times t₀ through t_n, and
the values y(1,1:2) through y(n+1,1:2) correspond to
vectors (w_0,1,w_0,2) through (w_n,1,w_n,2)

Running the code (with Matlab or Octave) produces the following:

fig_6.9 in text

fig_6.11

As written, the program euler2.m performs an animation (under Matlab) of the pendulum.

The modifications needed to produce figure 6.11, consist of

• replacing "trapstep" with "eulerstep",
• placing a "%" at the beginning of the each of the five lines

    xbob = cos(y(k+1,1)-pi/2); ybob = sin(y(k+1,1)-pi/2);
    xrod = [0 xbob]; yrod = [0 ybob];
    set(rod,'xdata',xrod,'ydata',yrod)
    set(bob,'xdata',xbob,'ydata',ybob)
    drawnow; pause(h)
  

• adding a line

  plot(t,y(:,1),t,y(:,2))
  

  
  after the end of the for loop

fig_6.12

fig_6.13

fig_6.16.jpg

Summary of Chapter 6

Sample exam questions and solutions
(Needs password from class)

438-538_sp11_prj1.html
Project 1: Buckling of a Circular Ring by Shooting

u(x,0) = f(x) for all a ≤ x ≤ b,

u(a,t) = l(t) and u(b,t) = r(t) for all t ≥ 0 .

Forward difference (in time):

Like Euler's method.

fig_8.2aa

Note the problem when the time step k is too large.
The size of the time step is limited by a stability requirement.

Stability analysis (Von Neumann) uses exact formulas for the eigenvalues of the matrix A taking the solution one step forwards in time. Need abs. value of all eigenvalues < 1

The eigenvectors of this matrix A are identical with the eigenvectors of a simpler tridiagonal matrix T

       [ 1  -1   0   0   ...   0]
       [-1   1  -1   0   ...   0]
T   =  [ 0  -1   1  -1   ...   0]
         . 
       [ 0  ...   0  -1    1  -1]
       [ 0  ...   0   0   -1   1]

λ_j = 1 - 2 cos(πj/(m+1)), j=1,..,m
for T.

Stability analysis leads to a constraint on σ = c k/h², where h is the step size in the space variable. To satisfy this constraint requires taking a small step size k in the time variable.

Like Backward Euler method.

fig_8.4aa.jpg

φ1 from top	φ2 from top
&phi1	&phi2

Sample files for Project 2:
gamip.txt
This file contains index pairs (i,j) for the network of web pages shown in Figure 12.1
If a pair (i,j) appearing in this file, then there is a link from page i to page j in the network.

gam.m Matlab/Octave program that reads the file gamip.txt
and writes a file gam.txt.

gam.txt
(Text) file containing adjacency matrix A.
A(i,j) = 1 if there is a link from page i to page j in the network,
A(i,j) = 0 if there is not such a link.
Each line of gam.txt corresponds to a row of A. The characters "1" and "0"
represent matrix elements and space " " separates elements.

gamnn.m
Matlab/Octave program that reads the file gam.txt
to determine the adjacency matrix A,
and then sums the elements in the rows of A
to construct a vector nn(i), i=1:n where nn(i) is the sum of row i of A.

hilberteig.m
Matlab/Octave program that sets up an example matrix H,
then computes the eigenvalues and eigenvectors of H,
determines the eigenvalue of largest absolute value, and
displays the corresponding eigenvector, where the eigenvector is normalized such that the sum of its elements is one.

Of interest (related to project 2 but not directly relevant)
The Second Eigenvalue of the Google Matrix
This theorem can be used to show that the condition |λ₂/λ₁| < 1 is satisfied,
(&lambda₁ = 1 by exercise 12.4.1 p. 562),
and so power iteration to find the first eigenvector p, should converge for most initial guesses.
(For each i, the component p_i is the probability of a web-surfer being on web page i.
The highest ranked pages will have the largest values of p_i.)

http://nptel.iitm.ac.in/courses/Webcourse-contents/IIT-KANPUR/Numerical%20Analysis/numerical-analysis/kadalbajoo/lec1/fnode3.html
A proof that the Maximum Absolute Column Sum Norm is the same as
the matrix norm ||A||₁ induced by the vector norm ||x||₁
(one way of working exercise 12.4.1b)