Download the files
num2bin.m
and bin2num.m (note that you
have to put the files where
Matlab can find them)- [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 - or use
x = A\b
which does the above steps automatically (but without the option to save the result of the LU-decomposition).
Errors for Linear Systems
divdiff computes the divided differences
d1=f[x1,...,xn],
d2=f[x2,...,xn],...,dn=f[xn]
Algorithm d=divdiff(x,y)
d := y
For k=1 to n-1
For i=1 to n-k
di :=
(di+1-di)/(xi+k-xi)
evnewt evaluates the interpolating polynomial at the point t
Algorithm v=evnewt(d,x,t)
v := d1
For i=2 to n
v := v·(t-xi) + di
Equidistant nodes on the interval [a,b] are given by
h = (b-a)/(n-1); xj = a + (j-1)h for j = 1,...,nThen
|(x-x1)...(x-xn)| ≤ hn(n-1)!/4 for x in [a,b]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 the QR decomposition ("economy version"):
[Q,R] = qr(A,0)
b = Q'*y
c = R\b- Matlab shortcut (actually uses QR decomposition):
c = A\y
We want to find a zero of the function f in the interval [a,b] where f(a)f(b)<0:
xs = fzero('f',[a,b])
where the function f is given by the m-file f.m.
Example with @-function:
f=@(x)sin(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.