NEWS:

## Course information

• Introduction: Errors, Taylor Theorem
• Error propagation, forward error analysis, numerical stability
• Examples for unavoidable error, actual error   (m-file)
• How to compute 1-cos(x)   (m-file)
• IEEE 754 machine numbers and machine arithmetic
• Solving Linear Systems
• How to solve a linear system in Matlab:
• [L,U,p] = lu(A,'vector'); % Gaussian elimination choosing the pivot candidate with largest absolute value
y = L\b(p)                % forward substitution with permuted vector b
x = U\y                   % back substitution
• Matlab shortcut
x = A\b
This performs the three steps from above, but you can't save L, U, p. Don't use this if you solve several linear systems with the same matrix.
• Example for Gaussian Elimination with Pivoting
• Errors for Linear Systems
For the 1-norm condition number use
c1 = norm(A,1)*invnormest(L,U)
For the ∞-norm condition number use
ci = norm(A,Inf)*invnormest(U',L')
For the Cholesky decomposition use invnormest(C',C). The m-file invnormest_simple.m explains the basic algorithm for estimating ||A-1||
• How to solve sparse linear systems in Matlab:
A sparse matrix is a matrix where most elements are zero. In this case it is much more efficient to use the special sparse data structures in Matlab. All operations like *, \, lu have special efficient algorithms for matrices in sparse format. In particular the backslash command x=A\b uses clever row and column permutations to minimize "fill-in".
• initializing a sparse matrix: A sparse matrix can be described by giving a list of the nonzero elements, as vectors iv (contains i-values), jv (contains j-values), av (contains matrix entries). Then a sparse matrix is initialized by
A = sparse(iv,jv,av);
Example: The full matrix given by A = [0 7 0; 0 0 8; 9 0 0] can be defined as a sparse matrix by
A = sparse([1 2 3],[2 3 1],[7 8 9])
• Solving a linear system with a sparse matrix:
• [L,U,p,q] = lu(A,'vector'); % Gaussian elimination with row and column permutations
% (chosen to minimize "fill-in"): L*U gives the permuted matrix A(p,q)
y = L\b(p)                  % forward substitution with permuted vector b
x = U\y; x(q) = x           % back substitution, then undo permutation given by q
• or just use
x = A\b
which does the above steps automatically (but without the option to save the result of the LU-decomposition).
• Cholesky decomposition in Matlab: (only works if A is symmetric positive definite)
• C = chol(A), then C'*C gives A
• for sparse matrix: find row/column permutation vector p to minimize "fill-in":
[C,k,p] = chol(A,'vector'), then C'*C gives the permuted matrix A(p,p)
(ignore the return value k which is zero for a positive definite matrix)
• Gaussian elimination WITHOUT/WITH pivoting, Cholesky decomposition
• Interpolation (updated March 13)
• Try out polynomial interpolation!
• Divided difference algorithm for distinct nodes x1,...,xn and given function values y1,...,yn:
divdiff computes the divided differences d1=f[x1,...,xn], d2=f[x2,...,xn],...,dn=f[xn]
Matlab code: divdiff.m
evnewt evaluates the interpolating polynomial at the point t
Matlab code: evnewt.m
• Example: Polynomial interpolation with equidistant and Chebyshev nodes     (m-file)
• Equidistant nodes on the interval [a,b] are given by
h = (b-a)/(n-1);    xj = a + (j-1)h for j = 1,...,n
Matlab: x = linspace(a,b,n)
Chebyshev nodes on the interval [a,b] are given by
xj = (a+b)/2 + cos(π(j-1/2)/n)·(b-a)/2 for j = 1,...,n.
The node polynomial has the bound
|(x-x1)...(x-xn)| ≤ 2 ((b-a)/4)n for x in [a,b]
• Piecewise polynomial interpolation This covers piecewise linear interpolation, piecewise cubic Hermite interpolation, complete cubic spline, not-a-knot cubic spline
You can ignore the formulas for the piecewise cubic Hermite interpolation and for the cubic spline interpolation. You just need to know how to compute this in Matlab, see right below.
• Piecewise cubic interpolation in Matlab: The column vector `x` contains the x-coordinates of the nodes, the column vector `y` contains the function values at the nodes, the vector `yp` contains the derivatives at the nodes. The vector `xi` contains the points where we want to evaluate the interpolating function, the vector `yi` contains the corresponding values of the interpolating function.
• Cubic Hermite interpolation: Download `hermite.m` . Then use ```yi = hermite(x,y,yp,xi)```
• Complete cubic spline: points are given by column vectors x, y. Slopes at left and right endpoint are sl and sr.
`yi = spline(x,[sl;y;sr],xi)`
• Not-a-knot cubic spline:``` yi = spline(x,y,xi)```
Example for piecewise cubic interpolation
• historical origin of "splines": "used principally by ship architects ... splines are thin, tapering pieces of lancewood or red pine" in "Treatise on Mathematical Drawing instruments" from 1868
• "Curve fitting" for data with measurement errors: Linear least squares method
• Nonlinear equations
• Nonlinear system: Example for Newton method and fsolve     (m-file)
• Solving a nonlinear equation with Matlab:
We want to find a zero of the function f in the interval [a,b] where f(a) and f(b) have different signs:
`xs = fzero(@f,[a,b])`
Here the function f is given in the m-file f.m. If the function f is defined by an @-expression use `xs=fzero(f,[a,b])`.
Example with @-function: ` f=@(x)sin(x); xs=fzero(f,[3,4])`
• Solving a nonlinear system in Matlab using fsolve: (without using the Jacobian)
Write an m-file f.m containing the function f:
function y = f(x)
y = ...
Then use x = fsolve(@f,x0) where x0 is the initial guess.
Example: See the bottom of the first page of Nonlinear equations
You can specify additional options using optimset, e.g. the tolerance for function value and the tolerance for x:
opt = optimset('TolFun',1e-5,'TolX',1e-5);
x = fsolve(@f,x0,opt);
There are many other options (e.g., display each iteration, or use the Jacobian) , type doc fsolve for more information.
• Fixed Point Iteration and Contraction Mapping Theorem (proof of Newton-Kantorovich theorem)
• Nonlinear least squares problem:
We have a function y=f(x) where x is a vector of length n and y is a vector of length N, with N≥n. We want to find x such that ||f(x)||2 is minimal.
Simple example with n=2, N=3: Let f(x) = [x1-.4, x2-.8, x12+x22-1]T, use initial guess [0,0]T
• Gauss-Newton method: gaussnewton.m
• How to use lsqnonlin (without Jacobian):
f = @(x) [x(1)-.4; x(2)-.8 ; x(1)^2+x(2)^2-1];
[x,resnorm2] = lsqnonlin(f,[0;0]) % resnorm2 is ||f(x)||2

There are many other options (e.g., display each iteration, or use the Jacobian) , type doc lsqnonlin for more information.
• Nonlinear least squares problem: Example for Gauss-Newton method and lsqnonlin
• Numerical Integration (updated 05/03)
• Numerical integration with Matlab:   Q = integral(f,a,b)
Example: Find the integral of exp(-x^2) from 0 to 1:
f = @(x) exp(-x.^2);   Q = integral(f,0,1)
Note:
• Write your function f(x) so that it works for a vector x and returns a vector f(x) of function values f(x1),...,f(xn): use .* ./ .^ instead of * / ^ in the definition of f.
Example: for I = 01e-x^2dx use
f = @(x) exp(-x.^2); % works for vectors x
Q = integral(f,0,1)
• If your function f(x) only works for scalars x you MUST use the option 'ArrayValued','True':
f = @(x) exp(-x^2); % works only for scalars
Q = integral(f,0,1,'ArrayValued',true)
• You can specify tolerances for the relative and absolute error:
Q = integral(f,a,b,'RelTol',1e-13,'AbsTol',1e-13)
01f(x)dx is approximated by w(1)*f(x(1))+...+w(n)*f(x(n))
abf(x)dx is approximated by r*(w(1)*f(t(1))+...+w(n)*f(t(n))) where r=(b-a)/2; t=(a+b)/2+r*x
Method from class. But GaussLegendre.m is more accurate and efficient.
01(x+1)α(1-x)βf(x)dx is approximated by w(1)*f(x(1))+...+w(n)*f(x(n))
ab(x-a)α(b-x)βf(x)dx is approximated by r^(1+alpha+beta)*(w(1)*f(t(1))+...+w(n)*f(t(n))) where r=(b-a)/2; t=(a+b)/2+r*x
• Gauss-Laguerre quadrature with α>-1, b>0
0xαe-xf(x)dx is approximated by w(1)*f(x(1))+...+w(n)*f(x(n))
0xαe-bxf(x)dx is approximated by b^(-1-alpha)*(w(1)*f(t(1))+...+w(n)*f(t(n))) where t=x/b
Examples for Gauss-Legendre and Gauss-Jacobi quadrature

## Matlab programming

We will use Matlab to see how various algorithms work.

• Gentle introduction to Matlab (also explains how to use the Matlab interface, how to publish)
• We will NOT use SYMBOLIC Matlab commands (like sym, syms, ezplot) in this class. (We will, however, use vpa in the first assignment.)
• Matlab Primer gives a concise summary of the most important Matlab commands (but it was written for an older version of Matlab)

How to hand in Matlab results for homeworks:

• You have to write an m-file for each problem. (Typing commands at the >> prompt and printing this out is not allowed)
• You have to hand in all m-files together with the generated output and graphics.
The best way is to use the publish command in Matlab (if you call additional m-files you have to print them out separately).
• Only print out values which are asked for in the problem. (Make sure you have a semicolon at the end of lines.)
Label all numerical output and graphics: How to label output and graphics in Matlab
• Hand in additional pages (this can be handwritten) which show the "paper and pencil work", and answer all the questions asked in the problem.

How to use the publish command in Matlab: