Assignment 5, problem 2: Trajectory of satellite in the earthmoon system with initial position (x_{1},0) (red dot) and initial velocity (0,v_{2}) (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:0012:00.
 practice problems for final exam for quadrature, ODEs
solution
 Final Exam: Tuesday, December 18, 10:30am12: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 \(\Acy\_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^32, or separate mfile 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 CO_{2} 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, LUdecomposition,
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
Additional Course Material (Required Reading)
 Review of Taylor's theorem with Example
 Introduction: Errors, Taylor Theorem
 Error propagation, forward error analysis, numerical
stability
 How to compute 1cos(x) for x close to zero
(generated with publish from this mfile)
How to analyze error for y=1cos(x) theoretically
 Useful functions available in all programming languages:
expm1(x) = e^{x}1, log1p(x) = ln(1+x),
hypot(x,y) = (x^{2} + y^{2})^{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 1norm condition number use
c1 = norm(A,1)*invnormest(L,U)
For the ∞norm condition number use
ci = norm(A,Inf)*invnormest(U',L')
The mfile 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 "fillin".
 initializing a sparse matrix: A sparse matrix can be described by giving a list of the nonzero elements, as
vectors iv (contains ivalues), jv (contains jvalues), 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 "fillin")
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 LUdecomposition).
 Example: Solving a problem with 2500 unknowns on a 50 by 50 grid, mfile
 Interpolation
 Interactive demonstration: Interpolation with polynomial and cubic spline
 Divided difference algorithm for nodes x_{1},...,x_{n} (different from each other) and given function values y_{1},...,y_{n}:
Download divdiff.m, evnewt.m.
d = divdiff(x,y) computes the divided differences
d_{1}=f[x_{1},...,x_{n}],
d_{2}=f[x_{2},...,x_{n}],...,d_{n}=f[x_{n}]
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
(mfile)
 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 xcoordinates 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)
 Notaknot 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 mfile 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^3x1; xs = fzero(f,[1,2])
 Nonlinear equations and nonlinear least squares problem
 Nonlinear system: Example for Newton method and fsolve (mfile)
 Solving a nonlinear system in Matlab using fsolve:
(without using the Jacobian)
Write an mfile 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',1e5,'TolX',1e5);
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 GaussNewton method and lsqnonlin
 Numerical integration (aka quadrature)
 Example for adaptive 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(x_{1}),...,f(x_{n}): 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',1e13,'AbsTol',1e13)
 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',1e14,'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\) (GaussLegendre) and
\(\rho(x)=(xa)^\alpha (bx)^\beta\) (GaussJacobi).
We can also use the infinite interval \([0,\infty)\) and \(\rho(x)=x^\alpha e^{x}\) (GaussLaguerre).
Matlab functions for finding \(x_1,\ldots,x_n\) and \(w_1,\ldots,w_n\):
 GaussLegendre 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=(ba)/2; t=(a+b)/2+r*x
 GaussJacobi quadrature
with \(\alpha,\beta>1\) : [x,w] = GaussJacobi(n,alpha,beta)
\(\displaystyle { \int_{1}^1 (x+1)^\alpha (1x)^\beta f(x)\,dx } \)
is approximated by w(1)*f(x(1))+...+w(n)*f(x(n))
\(\displaystyle \int_a^b (xa)^\alpha (bx)^\beta f(x)\,dx\)
is approximated by r^(1+alpha+beta)*(w(1)*f(t(1))+...+w(n)*f(t(n))) where
r=(ba)/2; t=(a+b)/2+r*x
 GaussLaguerre 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^(1alpha)*(w(1)*f(t(1))+...+w(n)*f(t(n))) where t=x/b
Examples for GaussLegendre and GaussJacobi quadrature
 Ordinary Differential Equations
 Matlab mfiles 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 tvalues, and ys is an array where the
jth column contains the values of the function y_{j}(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 y_{1}(T) = y(T)
plot(ts,ys(:,1)) % first colum of ys contains values of
y_{1}(t) = y(t)
Specifying desired accuracy:
opt = odeset('RelTol',1e5,'AbsTol',1e8);
[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.
 Matlab can be downloaded for free from TERPware.
 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 mfile for each problem. (Typing commands at the >> prompt and printing this out is not allowed)
 You have to hand in all mfiles together with the generated output and graphics.
The best way is to use the publish command in Matlab
(if you call additional mfiles 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.