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). - For large sparse matrices (most elements are zero):
use sparse to initialize matrix, e.g., A=sparse(n,n).
All operations like *, \, lu also work
with matrices in sparse format. In particular the backslash command
x=A\b uses clever row and column permutations to minimize "fill-in".
The most efficient way to initialize a sparse matrix is using a list of nonzero elements, given by vectors iv (contains i-values), jv (contains j-values), av (contains matrix entries):A = sparse(iv,jv,av);
E.g., the full matrix given by A = [0 7 0; 0 0 8; 9 0 0] can be defined as a sparse matrix byA = sparse([1 2 3],[2 3 1],[7 8 9])
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,...,nIn Matlab you can use x = linspace(a,b,n) to define a vector of equidistant nodes.
The node polynomial has the bound
|(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(The node polynomial has the bound(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)
In all three cases these we can first compute the spline as pp (a special structure for representing piecewise polynomials), and later evaluate the spline: E.g., for the not-a-knot spline:
pp = spline(x,y);
yi = ppval(pp,xi);
|| Ax - b || = minimal.
- 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])- Without using the Jacobian:
Write an m-file f.m containing the function f:function y = f(x)
Then use x = fsolve(@f,x0) where x0 is the initial guess.
y = ... - With using the Jacobian:
Write an m-file f.m containing the function and the Jacobian:function [y,A] = f(x)
Then use
y = ...
A = ...opt = optimset('Jacobian','on'); x = fsolve(@f,x0,opt);
where x0 is the initial guess. - You can specify additional options using optimset, e.g. the tolerance for the size of the
function value and the tolerance for the change in x:
opt = optimset('TolFun',1e-5,'TolX',1e-5);
There are many other options, type doc fsolve for more information.
x = fsolve(@f,x0,opt);
Q = w1f(x1) + ... + wnf(xn)with nodes x1,...,xn and weights w1,...,wn.
Download gauleg.m.
[x,w] = gauleg(n) gives n Gauss nodes x(1),...,x(n) and weights w(1),...,w(n) for the interval [-1,1]. For integration on an interval [a,b] use the nodes xt and weights wt given by
xt = (a+b)/2 + x*(b-a)/2; wt = (b-a)/2*w;
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,...). Note: The function f has to accept a vector x and return a vector y=f(x) of function values. Therefore use .* ./ .^ instead of * / ^ in the definition of f.
Example: y''+2y'+3y = sin(t), y(0)=2, y'(0)=-1; find y(t) for t in [0,5].
f = @(t,y) [y(2);-2*y(2)-3*y(1)+sin(t)];
[ts,ys] = Euler(f,[0,5],[2;-1],20);
plot(ts,ys(:,1),'.-'); hold all
[ts,ys] = Euler(f,[0,5],[2;-1],100);
plot(ts,ys(:,1),'.-'); hold off; legend('N=20','N=100')