Download the files
num2bin.m
and bin2num.m (note that you
have to put the files where
Matlab can find them) x = A\b or[L1,U] = lu(A); % compute
factorization A=L1 U using Gaussian elimination with pivoting,
choosing pivot candidates with largest absolute valuey = L1\b; % "forward substitution"
x = U\y; % back susbstitution
condition = condestdec(L1,U) %
estimate condition number of A (for 1-norm)norm(x,1) ,
norm(x,2) , norm(x,inf) condition numbers in Matlab:
cond(A,1),
cond(A,2), cond(A,3) (in all cases
inv(A) is computed)estimate for cond1(A):
condest(A)
(actually computes LU decomposition of A) If decomposition A=L1 U is already computed, use
condestdec(L1,U) (download
condestdec.m ,
for Matlab 5 use condestdec5.m
instead)The Chebyshev nodes on the interval [a,b] are given by
xj = (a+b)/2 + cos(Then(j-1/2)/n) (b-a)/2 for j = 1,...,n.
|(x-x1)...(x-xn)|2 ((b-a)/4)n for x in [a,b]
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 useyi = 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)
using normal equations:
- B = A'*A; z = A'*y; c = B\z
- using the QR decomposition ("economy version"):
[Q0,R0] = qr(A,0)
b0 = Q0'*y
c = R0\b0- Matlab shortcut (actually uses QR decomposition):
c = A\y
xs = fzero('f',[a,b]) for function in m-file
f.mExample with inline function:
f=inline('sin(x)','x'); xs=fzero(f,[3,4])[xg,wg] = findgauss(n) gives n Gauss nodes xg(1),...,xg(n) and weights wg(1),...,wg(n) . With tg = (a+b)/2 + xg*(b-a)/2 the integral of f(x) from a to b can be approximated by
Q = (b-a)/2*(wg(1)*f(tg(1)) + ... + wg(n)*f(tg(n))
The routine findgauss.m implements the method for finding the Gauss weights and nodes as explained in class. However, a more efficient and numerically stable way to compute this is given by gauleg.m.
[xg,wg] = gaussjac(n,alpha,beta) gives n Gauss nodes xg(1),...,xg(n) and weights wg(1),...,wg(n) . With tg = (a+b)/2 + xg*(b-a)/2 the integral of f(x)*(x-a)^alpha*(b-x)^beta from a to b can be approximated by
Q = ((b-a)/2)^(1+alpha+beta)*(wg(1)*f(tg(1)) + ... + wg(n)*f(tg(n)).
Example: gaussj_ex.m
f may be string f='x.^2', inline function f=inline('x.^2','x'), or anonymous function f=@(x) x.^2 . If the function is in an m-file f.m use quad('f',...) or quad(@f,...). Note: You MUST use .* ./ .^ instead of * / ^ in the definition of f.
ode45: Default
for RelTol is 1e-3, default for AbsTol is 1e-6
opts = odeset('RelTol',1e-4,'AbsTol',1e-7);
[t,y] = ode45('f',[t0,T],y0,opts);
ode*, error tolerance, event
locationode45: try this first (uses methods of order 4,5)ode23: sometimes better for lower accuracy (uses methods of
order 2,3)ode113: sometimes better for higher accuracy (uses variable
order) ode15s: for "stiff" problems, high accuracy (uses variable
order)ode23s: for "stiff" problems, low accuracy (uses order
2)