Assignment 5, problem 2: Trajectory of satellite in the earth-moon system with initial position (x1,0) (red dot) and initial velocity (0,v2) (blue arrow).
Drag the tip of the arrow with your finger/mouse and try to find periodic orbits!

NEWS:
• I will have office hours on Monday, Dec. 10, 10:00-12:00.
• practice problems for final exam for quadrature, ODEs
solution
• Final Exam: Tuesday, December 18, 10:30am-12:30pm in the usual room (MTH0411).
You can bring a "cheat sheet": letter size (8.5 by 11 inches), use front and back, MUST BE HANDWRITTEN.
The exam covers
• the topics of exam 1, exam 2 (see below)
• Ordinary differential equations: stable and unstable ODEs
Euler method, Improved Euler method
behavior of the errors: (i)  at fixed $$T$$ for $$h\to 0$$, (ii)  for fixed $$h$$ as $$t\to\infty$$
systems of ODEs, 2nd order ODEs
Matlab: [ts,ys] = ode45(f,[t0,T],y0)
• Assignment 5, due Monday, Dec. 10 at 9pm.
Note: I will drop the lowest of the 5 assignment scores.
• Exam 2 was on Tue, Nov. 27.
Exam 2 covered
• Linear least squares problem: minimize $$\|Ac-y\|_2$$, curve fitting, normal equations, Matlab shortcut c=A\y
• Nonlinear equation $$f(x)=0$$ for $$x\in\mathbb{R}$$: bisection method, secant method, Newton method, Matlab xs = fzero(f,[a,b])
• Nonlinear system $$f(x)=0$$ for $$x\in\mathbb{R}^n$$: Newton method, Matlab: xs=fsolve(f,x0)
• Numerical integration: midpoint, trapezoid, Simpson rule, composite rules, error formulas (I will print the error formulas for the basic rules on the exam), Matlab Q=integral(f,a,b) (if f returns multiple function values for vector x), Q=integral(f,a,b,'ArrayValued',true) works in general
• How to use functions in Matlab: either @-functions like f=@(x) x^3-2, or separate m-file f.m starting with function y=f(x)
• NOT ON EXAM: nonlinear least squares, adaptive quadrature, Gaussian quadrature
Practice problems for Exam 2 ,   solution
• Assignment 4, due Friday, Nov. 16 at 9pm, posted Nov. 6.
Hints
• Assignment 3, due Tuesday, Nov. 6 at 9pm, posted Oct. 26.
Hints
• How to import CO2 data from NOAA into Matlab
• Practice problems for Exam 1.
solution
• Exam 1 was on Tuesday, October 23.
Exam 1 covered
• errors and error propagation: machine epsilon, condition $$c_f(x)$$ of function, relative error for $$z=x+y$$, subtractive cancellation
• linear systems: Gaussian elimination without/with pivoting, LU-decomposition, residual vector, norms, errors, condition $$\mathrm{cond}(A)$$ of matrix
• interpolation with polynomials: divided difference method, error formula, case of multiple nodes
You need to know the Matlab commands \, lu, norm.
The following topics will NOT be on the exam
• interpolation with piecewise polynomials, splines
• least squares problem
• nonlinear equations
• Assignment 2, due Thursday, October 18 at 9pm, posted Oct. 9.
Hints
• Assignment 1, due September 6, posted Aug. 30.
Hints

Course information

• Review of Taylor's theorem with Example
• Introduction: Errors, Taylor Theorem
• Error propagation, forward error analysis, numerical stability
• How to compute 1-cos(x) for x close to zero     (generated with publish from this m-file)
How to analyze error for y=1-cos(x) theoretically
• Useful functions available in all programming languages:
expm1(x) = ex-1,     log1p(x) = ln(1+x),     hypot(x,y) = (x2 + y2)1/2
• Machine numbers and machine arithmetic
• 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.
• Errors for Linear Systems
• Estimating the condition number if L,U are known: Download invnormest.m.
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')
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")
y = L\b(p)                  % forward substitution with permuted vector b
x = U\y; x(q) = x           % back substitution, then undo permutation given by q
• Matlab shortcut
x = A\b
This performs the above steps automatically (but without the option to save the result of the LU-decomposition).
• Example: Solving a problem with 2500 unknowns on a 50 by 50 grid,     m-file
• Interpolation
• Interactive demonstration: Interpolation with polynomial and cubic spline
• Divided difference algorithm for nodes x1,...,xn (different from each other) and given function values y1,...,yn:
d = divdiff(x,y) computes the divided differences d1=f[x1,...,xn], d2=f[x2,...,xn],...,dn=f[xn]
p = evnewt(d,x,t) evaluates the interpolating polynomial at the point t (for a vector t this returns a vector p of values)
Example for interpolation
• Example: Polynomial interpolation with equidistant and Chebyshev nodes     (m-file)
• Try to find the optimal interpolation nodes: Interactive demonstration of the node polynomial
• Piecewise polynomial interpolation: p.w. linear, p.w. cubic Hermite, cubic spline
• 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
• Linear least squares problem
• 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 sign:
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]).
Examples with @-function:   f = @(x) sin(x); xs = fzero(f,[3,4])
f = @(x) x^3-x-1; xs = fzero(f,[1,2])
• Nonlinear equations and nonlinear least squares problem
• Nonlinear system: Example for Newton method and fsolve     (m-file)
• 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.
• Nonlinear least squares problem: Example for Gauss-Newton method and lsqnonlin
• Numerical integration (aka quadrature)
• Numerical integration with Matlab:   Q = integral(f,a,b)
Example: Find the integral of $$e^{(-x^2)}$$ from 0 to 5:
f = @(x) exp(-x.^2);     % Use .* ./ .^ instead of * / ^
Q = integral(f,0,5)
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. E.g., for $$f(x)=\sin(x^2)/x$$ use f=@(x)sin(x.^2)./x
• If your function f(x) only works for scalars x you MUST use the option 'ArrayValued',true:
Example: f = @(x) exp(-x^2); Q = integral(f,0,5,'ArrayValued',true)
• You can specify tolerances for the relative and absolute error:
Q = integral(f,a,b,'RelTol',1e-13,'AbsTol',1e-13)
• This also works for "improper integrals", e.g., $$\displaystyle \int_0^\infty x^{-1/3} e^{-x}\cos x\,dx$$
f = @(x) x.^(-1/3).*exp(-x).*cos(x)
Q = integral(f,0,Inf)
Q = integral(f,0,Inf,'RelTol',1e-14,'AbsTol',0)   % want relative error$$\le 10^{-14}$$
• Gauss quadrature: section 5 of Numerical integration gives examples for Gauss quadrature with several weight functions $$\rho(x)$$.
The most important cases for an interval $$[a,b]$$ are $$\rho(x)=1$$ (Gauss-Legendre) and $$\rho(x)=(x-a)^\alpha (b-x)^\beta$$ (Gauss-Jacobi).
We can also use the infinite interval $$[0,\infty)$$ and $$\rho(x)=x^\alpha e^{-x}$$ (Gauss-Laguerre).
Matlab functions for finding $$x_1,\ldots,x_n$$ and $$w_1,\ldots,w_n$$:
• Gauss-Legendre quadrature: [x,w] = GaussLegendre(n)
$$\displaystyle { \int_{-1}^1 f(x)\,dx }$$ is approximated by w(1)*f(x(1))+...+w(n)*f(x(n))
$$\displaystyle \int_a^b f(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
• Gauss-Jacobi quadrature with $$\alpha,\beta>-1$$ : [x,w] = GaussJacobi(n,alpha,beta)
$$\displaystyle { \int_{-1}^1 (x+1)^\alpha (1-x)^\beta f(x)\,dx }$$ is approximated by w(1)*f(x(1))+...+w(n)*f(x(n))
$$\displaystyle \int_a^b (x-a)^\alpha (b-x)^\beta 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 $$\alpha>-1,b>0$$ : [x,w] = GaussLaguerre(n,alpha)
$$\displaystyle { \int_0^\infty x^\alpha e^{-x}f(x)\,dx }$$ is approximated by w(1)*f(x(1))+...+w(n)*f(x(n))
$$\displaystyle \int_0^\infty x^\alpha e^{-bx} f(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
• Ordinary Differential Equations
• Matlab m-files for ODEs: dirfield.m (plot direction field), Euler.m (Euler method), IEuler.m (Improved Euler method)
• Example for differential equation:
Solving an initial value problem with ode45   ,   with Euler method   ,   with Improved Euler method
• Example for 2nd order differential equation:
Solving an initial value problem with ode45   ,   with Euler method   ,   with Improved Euler method
• Initial value problem for ODE: $\mathbf{y}' = \mathbf{f}(t,\mathbf{y}) , \qquad \mathbf{y}(t_0) = \mathbf{y}^{(0)}$ Solve with Matlab: [ts,ys] = ode45(f,[t0,T],y0);
Here ts is a column vector of t-values, and ys is an array where the j-th column contains the values of the function yj(t) for each time step.
Example: $$y''+y'+2y=t$$ with initial conditions $$y(0)=1, y'(0)=-2$$.
Print out $$y(10)$$ and plot $$y(t)$$ for $$t\in [0,10]$$.
• Rewrite as 1st order system: Let $$y_1(t):=y(t)$$, $$y_2(t):=y'(t)$$, then \begin{align} y_1' & = y_2 \\ y_2' & = t - y_2 -2y_1 \end{align}
• Matlab code:
f = @(t,y) [y(2); t - y(2) - 2*y(1)];
[ts,ys] = ode45(f,[0,10],[1;-2]);
ys(end,1)         % last row, 1st column of ys contains y1(T) = y(T)
plot(ts,ys(:,1))  % first colum of ys contains values of y1(t) = y(t)
Specifying desired accuracy:
opt = odeset('RelTol',1e-5,'AbsTol',1e-8);
[ts,ys] = ode45(f,[t0,T],y0,opt);

Solving a "stiff" problem: use ode15s in place of ode45.
Example: Solve the ODE $$y' = -10^5 y + \sin t$$ with initial condition $$y(0)=1$$ for $$t\in [0,10]$$.
f = @(t,y) -1e5*y + sin(t);
opt = odeset('Stats','on');       % display information about work
[ts,ys] = ode45(f,[0,10],1,opt);  % uses 301318 steps, 1.9 million evaluations of f
[ts,ys] = ode15s(f,[0,10],1,opt); % uses 96 steps, 201 evaluations of f

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, subs, diff, taylor etc.) in this class. (We will, however, use vpa in the first assignments to compare machine arithmetic with higher precision arithmetic.)
• Matlab Primer gives a concise summary of the most important Matlab commands. It was written for an older version of Matlab, so skip sections 14, 15, 17.

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 display and label numerical 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:
Use lines starting with %% for the main title and for each section title.
Lines starting with % right after a section title will be printed at the beginning of each section. Use this for your observations about the problem.