CSc 2262 Numerical Methods Fall 2004
More Notes
Fourier Series – Another View
f(x) over [-¥,¥] F(f) = ò-¥¥f(t)e-j2Pftdt with j=Ö-1
Note that f(t) = ò-¥¥F(f)ej2Pftdf and ejq=cos(q)+j*sin(q)
Discrete: x0, x1,…, xN-1 with xi – xreal+j* ximag
X(n) = (1/n)*Sk=1N-1
x(k)*e-jk2Pn/N Note
that x(n) = Sk=1N-1
X(k)*ejk2Pn/N for n=0,…,N-1
Magnitude = ||X(n)|| = (xreal*xreal + ximag*ximag)0.5
Phase = tan-1(ximag/xreal)
Properties
Linear – af(t)+bg(t) à aF(f)+bG(f)
Scaling – f(t/a) à aF(a*f)
Shifting – f(t+a) à F(f)e-2jPaf
Modulation – f(t)e-2jPaf à f(t-a)
Duality – Xk à (1/N)xN-k
Algorithms – Discrete, Fast Fourier
Transform
Examples
f(x)=1 if xe[0,P]; -1 if xe[P,2P) à
F(x) = (4/P)*[sin(x)+sin(3x)/3+sin(5x)/5+…]
f(x) = x if xe[-P,P] à F(x)=2*[sin(x)-sin(2x)/2+sin(3x)/3 …]
Example:
f(x) = 0 if xe(-P,0);P if xe[0,P] à
F(x) = (P/2)+2*[sin(x)+sin(3x)/3+sin(5x)/5+…]
Numerical Differentiation
Function f(x) continuous and differentiable over a given interval [a,b]
f ’(x) = derivative of f(x) = df(x)/dx
f ’(x) = limhà0 [f(x+h)-f(x)]/h
f ’(x) ~ [f(x+h)-f(x)]/h = Dhf(x) numerical derivative for small h (stepsize)
forward difference
Taylor series à f(x+h) = f(x) + hf ’(x) + (h2/2)*f ‘’(c) for cÎ[x,x+h]
à Dhf(x) = (1/h)*[f(x)+hf ’(x)+(h2/2)*f’ ’’(c)-f(x)] = f ’(x)+(h/2)*f ‘’(c)
à error = f ’(x) - Dhf(x) = -(h/2)f ‘’(c)
Example: f(x) = cos(x) Note that f ’(x) = -sin(x) and f ‘’(x) = -cos(x)
Let x = P/6 = 30° sin(x)=0.5000 and cos(x)=0.8660
Dhf(x) = [cos(x+h)-cos(h)]/h
Error = -(h/2)cos(c) à
error proportional to h but can’t let h get too small – why?
Backward Difference f ’(x) ~ [f(x)-f(x-h)]/h with hñ0 so error = (h/2)f ‘’(c)
Interpolating
Polynomial (Lagrange)
N = 2, x0 = x1-h, x1, x2 = x1+h yi=f(xi) for i=0,1,2
P2(x)=y0*(x-x1)(x-x2)/(x0-x1)(x0-x2)+y1*(x-x0)(x-x2)/(x1-x0)(x1-x2)+
y2*(x-x0)(x-x1)/(x2-x0)(x2-x1)
= y0*(x-x1)(x-x2)/2h2+y1*(x-x0)(x-x2)/(-h2)+y0*(x-x0)(x-x1)/(2h2)
P2’(x) = (2x-x1-x2)y0/2h2+(2x-x0-x2)y1/(-h2)+(2x-x0-x1)y2/2h2
à P2’(x1) = (x1-x2)y0/2h2+(2x1-x0-x2)y1/(-h2)+(x1-x0)y2/2h2 = [f(x2)-f(x0)]/2h
à f’(x1) ~ P2’(x) = [f(x1+h)-f(x1-h)]/2h = Dhf(x1) Central difference formula
f’(x)-Pn’(x) = Yn(x)*f(n+2)(c1)/(n+2)! + Yn’(x)*f(n+1)(c2)/(n+2)! With Yn(x)=Pj=1n(x-xj)
c1 and c2Î[min xj value, max xj value] and the xj values include x for this interval
n=2 yields Yn(x)=(x-x0)(x-x1)(x-x2) and
Yn’(x)=(x-x1)(x-x2)+(x-x0)(x-x2)+(x-x0) (x-x1)
à Yn’(x1)=(x1-x0)(x1-x2) = -h2
à f’(x1) ~ [f(x1+h)-f(x1-h)]/2h = (h2/6)*f”’(c2)
Undetermined
Coefficients
f ‘’(x)~Dh(2)f(x)=Af(x+h)+Bf(x)+Cf(x-h)
Taylor à f(x-h)~f(x)-hf ’(x)+h2/2*f ‘’(x)-h3/6*f ‘’’(x)+h4/24*f ‘’’’(x)
à f(x+h)~f(x)+hf ’(x)+h2/2*f ‘’(x)+h3/6*f ‘’’(x)+h4/24*f ‘’’’(x)
à Dh(2)f(x) ~ (A+B+C)*f(x) + h*(A-C)*f ’(x) h2/2*(A+C)*f ‘’(x) +
h3/6*(A-C)*f ‘’’(x)+h4/24*(A+C)*f ‘’’’(x)
Dh(2)f(x) ~ f ”(x) à (A+B+C) = coefficient of f(x) = 0;
h(A-C) = coefficient of f ’(x) = 0;
h2/2*(A+C) = coefficient of f ‘’(x) = 1
à A = C = 1/h2 and B = -2/h2
à Dh(2)f(x)
~ [f(x+h)-2f(x)+f(x-h)]/h2
à f ‘’(x)- Dh(2)f(x) ~ -h2/12*f ‘’’’(x)
Newton-Cotes Formulae
I(f) = òabf(x)dx ~ òabPn(x)dx
where Pn(x) = åj=0najxj = a polynomial
P1(x) = [(b-x)f(a)+(x-a)f(b)]/(b-a) -->
T1(f) = òab[f(a)+{f(b)-f(a)}*(x-a)/(b-a)]dx = (b-a)[f(b)+f(a)]/2
Error = (1/12)f ‘’(c)(b-a)3
Example: f(x) = 1/(1+x) a=0 and b=1 Note: F(x) = ln(1+x)
T1 = 0.75 but I(f) = ln(2) = 0.693147
Composite: n intervals h=(b-a)/n xj=a+h*j j=0,1,…,n
I(f) = òajb f(x)dx =åj=0n-1òxjxj+1 f(x)dx ~ Tn(f) = (h/2)åj=0n-1[f(xj)+f(xj+1)]
= (h/2)[f(x0)+f(xn)+2åj=1n-1f(xj)]=(b-a)*[f(x0)+f(xn)+2åj=1n-1f(xj)]/2n
where (b-a) = width and [f(x0)+f(xn)+2åj=1n-1f(xj)]/2n = average heigth
Error -(b-a)3åj=1n f ‘’(cj)/12n3
Example: n = 4
h=(b-a)/4 x0=a, x1=a+h=(3a+b)/4, x2=a+2h=(a+b)/2, x3=a+3h=(a+3b)/4, x4=a+4h=b
f(x) = e^(-x2), n=4, a=0, b=1 --> T4=0.7468241…
f(x) = 1/(1+x2), n=4, a=0, b=4 --> T4=1.32581766…
I = tan-1(4)
f(x) = 1/(2+cos(x)), n=4, a=0, b=2p --> T4=0.7468241…
I=2p/Ö3
Simpson’s Rule
N=2 c=(a+b)/2 h=(b-a)/n
I(f) ~ òabP2(x)dx
= òab
[(x-c)(x-b)f)a)/(a-c)(a-b)+
(x-a)(x-b)f(c)/(a-c)(b-c)+(x-a)(x-c)f(b)/(b-a)(b-c)]dx
=òab
[(x-c)(x-b)f)a)/(a-c)(a-b)]dx+òab[(x-a)(x-b)f(c)/(a-c)(b-c)]dx+
òab
[(x-a)(x-c)f(b)/(b-a)(b-c)]dx
Note:
òab(x-c)(x-b)/(a-c)(a-b)dx
= (1/2h2)òaa+nh(x-c)(x-b)dx
=(1/2h2)òa2h(u-h)(u-2h)du with
u=x-a
=(1/2h2)[u3/3-3u2h/2+2h2u]02h
= h/3
à I(f) ~ S2(f) = (h/3)[f(a)+4f(c)+f(b)]
Example: a=0 b=1 f(x) = 1/(1+x) c=(a+b)/2 = 1/2 h=(b-a)/2 = 1/2
S2(f) = [(1/2)/3] [(1+4(2/3)+1/2] = 25/36 = 0.694444
I(f) = ln(2) = 0.693147
Simpson’s rule works well if f(x) nearly quadratic
Composite: n intervals h=(b-a)/n xj=a+h*j j=0,1,…,n
I(f)=òabf(x)dx=òx0xn f(x)dx=òx0x2f(x)dx +òx2x4f(x)dx +…+òxn-2xn f(x)dx
Sn(f) =
(h/3)[f(x0)+4f(x1)+2f(x2)+4f(x3)+2f(x4)+…+2f(xn-2)+4f(xn-1)+f(xn)]
=(b-a)[f(x0)+f(xn)+4åj=1,j oddn-1f(xj)+2åj=1,j evenn-1f(xj)]/6
Error = [(b-a)5 /180n4]f(4)
where f(4 = average f(4) value
Composite’ use
of 3rd order Lagrange polynomial
I(f) ~ (3h/8) [f(x0)+3f(x1)+3f(x2)+
f(x3]
=(b-a) [f(x0)+3f(x1)+3f(x2)+
f(x3]/8
Error
= -[3h5/80]*f(4)(c) = -[(b-a)5/6480]*f(4)(c)
Can do higher orders
Gaussian Quadrature
Assume a = -1 and b = 1 I(f) = òabf(x)dx = ò-11f(x)dx ~ In(f) = Sj=1nwjf(xj)
Choose x and w values so In(f) = I(f) for as large a degree n as possible for simple functions
n = 1
f(x) =1 à I(f) = ò-11dx = 2 = I1(f) = w1 à w1 = 2
f(x) = x à I(f) = ò-11xdx = 0 = I1(f) = w1x1 = 2x1 à x1 = 0
à I1(f) = 2f(0)
Note: I(f) = òabf(x)dx ~ òabf[(a+b)/2]dx = (b-a)f[(b+a)/2]
a = -1 and b = 1 à I(f) ~ 2f(0) midpoint formula
n = 2
I2(f) = w1f(x1)+w2f(x2)
f(x)=1,x,x2,x3 à w1+w2=2; w1x1+w2x2=w1x13+w2x23=0; w1x12+w2x22=2/3
à w1=w2=1, x1=-Ö3/3, x2=Ö3/3 à I2(f) = f(-Ö3/3)+f(Ö3/3)
n>2
In(f) = Sj=1nwjf(xj)
f(x)= 1,x,x2,…,x2n-1 à Sj=1n wj=0=Sj=1n wjfxj3, 2/3=Sj=1n wjxj2, … ,
Sj=1n wjxj2n-1=2/(2n-1)
nonlinear equations
If a¹-1 and/or b¹1 then let x = [b+a+(b-a)t]/2 where tÎ[-1.1]
à I(f)= (b-a)/2ò-11f([b+a+(b-a)t]/2)dt
Weighted Gaussian
Quadrature
Replace x with w(x); e.g., 1/Ö(1-x2) or log(1/x) or Öx or 1/Öx
Let w(x) = 1/Öx and a=0 and b=1
n = 1
f(x) =1 à I(f) = ò01(1/Öx)dx = 2 = I1(f) = w1 à w1 = 2
f(x) = x à I(f) = ò01(x/Öx)dx = 2/3 = I1(f) = w1x1 = 2x1 à x1 = 1/3
à I1(f) = 2f(1/3)
n = 2
I2(f) = w1f(x1)+w2f(x2)
f(x)=1,x,x2,x3 à w1+w2=2; w1x1+w2x2=2/3; w1x12+w2x22=2/5; w1x13+w2x23=2/7
à w1=1+Ö30/18; w2=1-Ö30/18; x1=3/7-2Ö30/35; x2=3/7+2Ö30/25
n>2
Sj=1n wj=2; Sj=1n wjfxj=2/3; Sj=1n wjxj2=2/5; … ; Sj=1n wjxj2n-1=2/(4n-1)
Word problem: A farmer has some ducks, all of which weigh the same, and some ducklings, all of which weigh the same. He finds that 3 ducks and 2 ducklings together weigh 32 kg.(kilograms), while 4 ducks and 3 ducklings together weigh 44 kg. Help the farmer figure out how much 2 ducks and one duckling weigh together.
Let x=weight of a duck and y=weight of a duckling
3x+2y=32 and 4x+3y=44 à x=8 and y=4
à 2x+y=20
Equations: In general, we have åj=1naijxj=bi for i,…m (m equations) and j=1,…n (n unknowns – x) – Note: the a’s and b’s are real constants
Note: if m<n, too few equations à no unique solution (later, we shall see that m may be reduced if dependencies exist)
Example: Let aij = max(i,j), bi=1 for all i,j à xj =0 for j=1,..,n-1 and xn=1/n. If n = 3, we have a11=1, a12=2, a13=3, a12=a22=2, a23=3, a13=a23=a33=3, b1=b2=b3=1 à x1=x2=0 and x3=1/3.
Linear Algebra - Matrix arithmetic
Scalar – number
Vector – one-dimensional array of numbers
Matrix – two-dimensional array of numbers
Matrix AmXn
Transpose
BnXm=AnXmT à bij= aji
for all i,j
Scalar multiplication - BmXn=sAmXn à bij= s*aij for all i and j
Addition (Subtraction) - CmXn=AmXn±BmXn à cij=aij±bij for all i and j
Multiplication - CmXs=AmXn±BrXs à cij=åk=1naik*bkj for all i and j and where n = r
Note: cij = dot or vector product of ith row of A (a vector) and jth column of B (a vector)
Back to equations:
Let A = matrix with m rows and n columns = {aij}mCn, b = (column) vector with m rows = {bi}mC1, and x = (column) vector with n rows = {xj}nC1. This means that our equation set reduces to Ax=b. For our example problem, A=[1 2 3; 2 2 3; 3 3 3] and b=[1;1;1] à x=[0;0;1/3]
Assume that m = n. Now, n is called the order of the equation set or of the matrix A or of the vectors b and x. If m=n, the matrix A is called square. The elements {aii} are called the main diagonal. If a square matrix has ones for the main diagonal and zeros elsewhere, it is called I, the identity matrix. Obviously, Iij=δij (the Kronecker delta). Note that A*I = I*A = A. We shall see that sometimes a matrix A-1 or A inverse, exists, such that A-1*A = A* A-1 = I. The purpose of this will be seen as if A-1 exists, and we have A*x=b, we can find x via A-1*A*x = I*x = x = A-1*b. Note that the determinant det(A)¹0 (A is nonsingular) for A-1 to exist. Moreover, the inverse A-1 is unique. For what it is worth, there is the zero matrix, 0, which has all zero elements.
Rules
(A+B)+C=A+(B+C); (AB)C=A(BC); A+B=B+A; A(B+C)=AB+AC; (A+B)C=AC+BC; (AB)T=BTAT; (A+B)T=AT+BT; (cA)-1=(1/c)A-1; (AB)-1=B-1A-1; det(AB)=det(A)det(B); det(AT)=det(A); det(cA)=cndet(A) with n=order(A) and c=constant
Assume that n>0 and A is a square matrix of order n with det(A)¹0. Then for each value of vector b, there is a unique solution x to Ax=b; note that this implies that there is at least one solution. Moreover, AB does not always equal BA (if there is equality, the matrices are said to commute). And if b=0, x=0.
Elementary row operations
How did we solve the duck/duckling word problem? Consider A=[3 2; 4 3] and b=[32; 44]. A-1 = [3 –2; -4 3]; x=A-1*b=[8; 4]. Note that A-1*A=I and that det(A) = a11a22- a21a12a=3*3-4*2=9-8=1¹0. To get the solution, we need elementary row operations.
i) interchange two rows (e.g., row 2 put
above row 1);
ii) multiply a row by a nonzero scalar (e.g.,
multiply row 2 by –3);
iii) add a nonzero multiple of one row to
another row (e.g., add 3*row 1 to row 3).
Don’s secret dirty trick
If you do an elementary row operation to I of same size as A to get I’
and multiply I’*A, you get A’ which is A with the same elementary row operation
done to it.
Example: A = [1 2 3; 2 2 3; 3 3 3] and n=3.
To interchange rows 1 and 2, do this to I =[1 0 0; 0 1 0; 0 0 1] to get
I ’= [0 1 0; 1 0 0; 0 0 1] and you can see that I ‘A=A ’=[2 2 3; 1 2 3; 3
3 3].
To now multiply row 2 of A ’ by –3, do this to I to get I ’’=[1 0 0; 0 –3 0; 0 0 1] and find A ’’= I ‘‘*A’ = [2 2 3; -3 –6 –9; 3 3 3].
To now add 3 * row 1 to row 3, do this this to I to get I ‘’’= [1 0 0; 0 1 0; 3 0 1] and find A ‘’’=I ‘’’*A ‘’ = [2 2 3; -3 –6 –9; 9 9 12].
Note that we could also do I ‘’’’= I ‘’’‘I ‘‘’I ‘= [0 1 0; -3 0 0; 0 3 1] and A ‘’’= I ‘’’’A
So, consider the extended matrix E=[A|b|I] where you take matrix A, add a column at the right and put b there, then add n columns to the right and put I there. Now do a series of elementary row operations on E until you get [I|x|A-1]. This is evident since you are multiplying E by A-1 to reduce A to A-1A=I, thus reducing b to A-1b = x, and A-1I to A-1. The trick is to determine the sequence of elementary row operations.
Gaussian elimination
Suppose n=3, A=[1 2 1; 2 2 3; -1 –3 0] and b=[0; 3; 2]. We can eliminate x1 from the
second and third equations by subtracting 2*equation 1 from equation 2 and by
subtracting –1*equation 1 from equation 3.
This yields A ‘=[1 2 1; 0 –2 1; 0 –1 1] and b’=[0; 3; 2]. Now, to eliminate x2 from the
third equation, subtract ½*equation 2 from equation 3 to yield A ‘’=[1 2
1; 0 –2 1; 0 0 ½] and b ‘’=[0; 3;
½]. The elimination steps force A ‘’ to
be upper triangular (consider the appearance) and now we can use substitution
to solve directly. After all, if
equation three is ½* x3=1/2, then x3=1. So, equation 2 gives us –2x2+ x3=-2x2+1=3
or x2=-1. This, in turn, can
be used in equation 1 to give x1+2x2+x3= x1+2(-1)+1=
x1-1=0 or x1=1.
Thus, x=[1; -1; 1]. Note that we
have used nothing but the elementary row operations to solve this.
To state this
algorithmically, consider original equation E(i) as åj=1naij(1)xj=bi(1)
for i=1,…n. Now, for k=1,2,…,n-1 do the
elimination steps: i) define multipliers mik=aik(k)/akk(k)
assuming akk(k)¹0 for i=k+1,…,n; ii)
subtract mik*E(k) from E(i) to eliminate xk from E(i), yielding aij(k+1)=
aij(k)-mikakj(k) for
i,j=k+1,…n and bi(k+1)=bi(k)-mik
bk(k) for I=k+1,…n.
Note that after these n-1 steps, the revised matrix A (i.e., the system
of equations) is in upper triangualar form.
Now, letting uij=aij(n) and gi=bi(n),
use back substitution to solve, giving i) xn= gn/unn;
and ii) xi=[gi-åj=i+1nuijxj]/uii for i=n-1,…,1.
One problem is what to do if akk(k)=0, violating the assumption and forcing division by zero (be careful that error does not force one into not realizing that the assumption is really violated but the variable is not exactly zero so one is falsely confident in the assumption). One can use partial pivoting, i.e., interchanging rows.
Basically, consider the E matrix where E=[A|b|I]. Note that E has m rows and 2n+1
columns. Note that eij = aij
for j=1,…,n; ein+1 = bi, and eijj=δij-n-1
(Kronecker delta) for j=n+2,,,,2n+1.
Now, the algorithm.
1.Consider a pivot row r and pivot column s, where r and s Î[1,2,…,n] Usually, we start with r=s=1.
2. If ers¹0 then we look
for eis¹0 for i¹r. If
no such i exisits, that means column s is all zeros, so det(A)=0 and no unique
solution exists. We can ignore this
column and continue the algorithm to see what happens if we wish or we can
terminate with an error message. If such an i exists, we interchange row r and
row i.
3. Now, we multiply row r by the 1 over the pivot element ers
so e ‘rj=(1/ers)*erj for j=1,…,2n+1. Note that this means that e’rs=1,
which is what we want.
4. Now, for each row i¹r we multiply the modified row r by -eis
and add it to row i, i.e., e’ij=(-eis)*e’rj+eij=(-eis)*(erj/ers)+eij
for j=1,2,…,2n+1. Note that eis=0
for i=1,…,r-1,r+1,…2n+1, which is what we want. Thus, when we are done with this iteration, we have column s as
an identity column, which is what we want.
5. Now we repeat steps 2-4 for r=s=2,3,…,n. Note that once you are done, A will become I, b will become x,
and I will become A-1. Thus
E=[A|b|I] à E ’= [I|x|A-1].
Example: ducks and ducklings).
Recall that a farmer has some ducks, all of which weigh the same, and some ducklings, all of which weigh the same; so he finds that 3 ducks and 2 ducklings together weigh 32 kg.(kilograms), while 4 ducks and 3 ducklings together weigh 44 kg; however, he needs help in figuring out how much 2 ducks and one duckling weigh together. Let x=weight of a duck and y=weight of a duckling à 3x+2y=32 and 4x+3y=44 à x=8 and y=4 à 2x+y=20.
A=[3 2; 4 3], b=[32; 44], x=[x1;x2], and I=[1 0;0 1]. Here n=2.
Note that det(A)=3*3-2*4=1¹0, so a unique solution exists.
E=[A|b|I]=[3 2 32 1 0; 4 3 44 0 1].
Note that E has n=2 rows and 2n+1=5 columns.
Let r=s=1 and we see that ers=e11=3¹0. Thus, E becomes [1 2/3 32/3 1/3 0; 4 3 44 0
1], obtained by multiplying row 1 by 1/ers. Now, multiplying the new row 1 by the
element in the pivot column but in the next row i=2, –eis=-e21=-4 and adding that to row I makes
the new E=[1 2/3 32/3 1/3 0; 0 1/3 4/3 –4/3 1]. Now, let r=s=2, so the new pivot element is ers=e22=1/3. Thus, multiplying row r=2 by 1/(1/3) lets E
become [1 2/3 32/3 1/3 0; 0 1 4 –4 3].
Finally, multiplying the new row 2 by –e12=-2/3 and adding to
row i=1 give the final E=[1 0 8 3 –2; 0 1 4 –4 3]=[I|x|A-1]. Note that x=[8;4] and A-1=[3 –2;
-4 3].
Note that the reverse sequence of modified identify matrices generated by
the elementary row operations is [1 –2/3; 0 1] [1 0; 0 3]; [1 0; -4 1] [1/3 0;
0 1] which (going right to left) multiplies row 1 by 1/3, multiplies row 1 by
–4 and adds to row 2, multiplies row 2 by 3, and multiplies row 2 by –2/3 and
adds to row 1. If we multiply these
four modified identity matrices together, we get A-1. Note that det(A-1)=1. Also, A-1A=A A-1=I and
A-1b=x=[8; 4].
LU Factorization
Ax=b à Ux=g where uij=aij(i) if j³i, 0 if j<i U upper triangular
U can be found using Don’s secret dirty trick but only downward (j<i) below main diagonal
Consider again Gaussian elimination where Sj=1naij(1)xj=bi(1) with mij=aij(j)/ajj(j) "i=j+1,…,n à aij(k+1)=aij(k)-mikakj(k) "i,j=k+1,…,n and bi(k+1i)=bi(k)-mikbk(k) "i=k+1,…,n
We can define L via lij=mij(i) if j<i, δij if j³i L lower triangular
Note that A=LU so this is called LU factorization
Example
n=4 A=[4 3 2 1; 3 4 3 2; 2 3 4 3; 1 2 3 4] and b[1; 1; -1; -1]
m21=3/4, m31=1/2, m41=1/4 so A ‘=[4 3 2 1; 0 7/4 3/2 5/4; 0 3/2 3 5/2; 0 5/4 5/2 15/4]. Now, m32=6/7 and m42=5/7 so A ‘’=[4 3 2 1; 0 7/4 3/2 5/4; 0 0 14/7 10/7; 0 0 10/7 20/7]. Now, m43=5/6 so A ‘’’= [4 3 2 1; 0 7/4 3/2 5/4; 0 0 12/7 10/7; 0 0 0 5/3] = U.
L = [1 0 0 0; ¾ 1 0 0; ½ 6/7 1 0; ¼ 5/7 5/6 1].
Note that L*U=A.
Matlab – lu function
x=A\b
Tridiagonal Matrix
Diagonal matrix - aij=0 if i¹j
Tridiagonal matrix - aij=0 if |i-j|>1
Let aii=bi, aii+1=ci, and aii-1=ai so uii=bi and uii+1=ci and lii=1 and lii-1=ai
à b1=b1, ajbj-1=aj, ajcj-1+bj=bj àb1=b1, aj=aj/bj-1, bj=bj–ajcj-1 "j=2,…,n
Consider Ax=f à Lg=f, Ux=g à g1=f1, gj=fj-ajgj-1 "j=2,3,…,n
à xn=gn/bn, xn-1j=[gj-cjxj+1]/bj "j=n-1,…,1
Iterative Methods Ax=b
Jacobi Iteration
Ax=b à Sj=1naijxj=Sj=1,j¹inaij(1)xj+aiixi=bi(1)bi à xi=(1/aii)[bi-Sj=1,j¹1naijxj] with aii¹0
Let x0=initial estimate for vector x. Then, xi(k+1)=(1/aii)[bi-Sj=1,j¹1naijxj(k)]
Gauss-Seidel
Similar to Jacobi but xi(k+1)=(1/aii)[bi-Sj=1,j>1naijxj(k)--Sj=1,j<1naijxj(k+1)]
*Advanced Stuff
Ax=b Nx=b+Px where A=N-P (splitting)
N nonsingular, can be diagonal or triangular or tridiagonal
elected to solve Nz=f easily
Nx(k+1)=b+Px(k) iteratively
Errors Stability Residual correction
Eigenvalues and Eigenvectors
A square matrix of order n
Av=lv à (lI-A)v=0
Obviously v=0 solves this but to get a solution v¹0 we need a
non-unique solution situation, i.e., A must be singular, i.e., f(l)=det(lI-A)=0.
Note that f(l) is a polynomial in the variable l, called the characteristic polynomial. The roots are called eigenvalues.
Consider a square matrix
A of order n.
det(A) = |A| =å j=1naijCij where Cij = a cofactor =
(-1)i+jMij
where Mij = a minor = det(A
without row i nor column j)
Examples: n=2 A={a b; c d} det(A)
= a(d)-b(c)
n=3
A = {a1 a2 a3; b1 b2 b3; c1 c2 c3}
det(A)=a1*det(b2 b3; c2 c3})–a2*det{b1 b3; c1
c3})+
a3*det({b1 b2; c1 c2})
=
a1*(b2c3-b3c2)-a2*(b1c3-c1b3)+a3*(b1c2-b2c1)
det(cA)=c for scalar c det(AB)=det(A)det(B)
det(A)det(A-1)=det(I)=1 à det(A-1)=1/det(A)
det(BAB-1)=det(B)det(A)det(B-1)=det(B)det(A)/det(B)=det(A)
det(B-1AB-cI)=det(B-1AB-B-1cIB)=det(B-1(A-cI)B))
=det(B-1)det(A-cI)det(B)=det(A-cI)
If A = U = upper
triangular matrix, then det(A)=P j=1najj
If 2 rows interchanged,
determinant changes sign
If row (or column) multiplied by c,
determinant multiplied by c
If any row (or column) = 0, determinant = 0
If any row (column) equals another row
(column), determinant = 0
If you multiply a row (column) by c and add
to another row (column),
determinant unchanged
Example: n=2. A={aij} and lI-A=[l-a11
-a12; -a21 l-a22] so f(l)=l2-(a11+a22)l+a11a22-a21a12. Let A=[1.25
0.75; 0.75 1.25]. Note det(A)=4. f(l)=l2-2.5l+1=0 à l1=0.5 and l2=2 as roots or eigenvalues.
Solving for vj,
the eigenvector (a value of v to
solve equation for a given eigenvalue) for eigenvalue lj, is difficult as eigenvectors are not
unique. For example, for nonzero
constant c, cv(j) is an eigenvector if v(j) is. However, we note that v(1)=[1; 1]
and v(2)=[-1;1] will work.
Example: n=3. A=[-7 13 -16; 13
-10 13; -16 13 -7]. f(l)=l3+24l2-405l+972=0 so l1=-36, l2=9, and l3=3. Suppose we select l=l1=-36.
Consider the original equation (lI-A)v=0. We have (lI-A)=[-29 -13 16; -13 -26 -13; 16 -13
-29]. To solve for the eigenvector,
suppose we let v1(1)=1 arbitrarily. Then, -13v2(1)+16v3(1)=29,
-26v2(1)-13v3(1)=13, and -13v2(1)-29v3(1)=
-16. Solving, we can get v2(1)=-1
and v3(1)=1.
For a symmetric matrix,
where AT=A or aij=aji for all i and j, the
eigenvalues are all real, the eigenvectors are mutually perpendicular (v(r
)*v(s)T=0 for r¹s), the length of the eigenvectors are 1 (||v(r)||=ÖSj=1n v(r) 2||=1). Moreover, for any vector x with n elements,
there are unique constants cj such that x=Sj=1n cjv(j) and ci=Sj=1n xjvj(i) .In addition, if we define U={v(j)}
then D=UTAU = {ljδij} and UUT=UTU=I. We sadly note that for nonsymmetric
matrices, we can get complex numbers for the eigenvalues and eigenvector
elements, so that the situation can get much more complicated.
Matlab
lambda=eig(A) gets
eigenvalues
[V,D]=eig(A) gets
eigenvectors as columns of V and eigenvectors as diagonal elements of D
Calculate numerically the
largest eigenvalue - Power method
Assume |l1| > |l2| ³…³|ln|.
This method calculates l1 and v(1). Let z(0)=initial
guess for v(1), usually done randomly. Define w(1)=Az(0) and a1=largest component of w(1). Now, iteratively, w(m)=Az(m-1),
am= maximum component of w(m), z(m)=w(m)/
am.
Moreover, l1(m)=wk(m)/zk(m-1)
where components k of both w and z are nonzero and usually the maximal
component for z(l) for some sufficiently large l. This converges to l1as mà¥.
Nonlinear equation systems
f(x,y)=0, g(x,y)=0
Newton’s method
(x0,y0)=initial
guess
r(x,y)=f(x0,y0)+(x-x0)fx(x0,y0)+(y-y0)fy(x0,y0)
where fx(x,y)=δf(x,y)/δx and fy(x,y)=δf(x,y)/δy. Note that r(x,y) = equation of tangent plane to z=f(x,y)
at (x0,y0,f(x0,y0)).
q(x,y)=g(x0,y0)+(x-x0)gx(x0,y0)+(y-y0)gy(x0,y0)
where gx(x,y)=δg(x,y)/δx and gy(x,y)=δg(x,y)/δy. Note that q(x,y) = equation of tangent plane
to z=g(x,y) at (x0,y0,g(x0,y0)).
Need to solve r(x,y)=0 and
q(x,y)=0 à x1=x0+δx and y1=y0+δy. Repeat iteratively.
Generalized method Modified
Method
Ordinary Differential
Equations(ODE)
Y’(x)=dy/dx=f(x,Y(x) x³(x0) where Y(x)=function of x
First order – first-order derivative
Practical Example
Newton’s law of cooling: Y ‘(x) = -k(Y(x)-A) where k is a positive constant, x is time, A is the temperature of the surrounding medium, and Y is the temperature of a cooling object – note that the rate of change of temperature is proportional to the difference in temperatures of the object and the surrounding medium.
Theory
Y ’(x)=g(x) à Y(x)=òg(x)dx+c c = arbitrary constant found with Y(x0)= Y0
Example: Y ’(x)=sin(x)à Y(x)= - cos(x)+c x0=P/3 and Y(x0)=2 à c=2.5
à Y(x)=2.5 - cos(x)
Y ’(x)=f(x,Y(x))
Y ‘(x) = a(x)Y(x)+b(x) with a, b continuous
f(x,z)=a(x)z+b(x)
Method of integrating factors
Y ‘(x)=lY(x)+b(x) x³ x0
à d(e-lxY(x))/dx=e-lxb(x) à e-xY(x)=c+òx0xe-ltb(t)dt
à c=e-lx0Y(x0)=e-lx0Y0
Example: Y ‘(x)=-[Y(x)]2+Y(x) Y0=0 à Y(x)=1/(1+ce-x) or Y(x)=0
If f(x,z) and df(x,z)/dz are continuous functions of x and z for all points in a neighborhood of (x0,Y0), then $ unique function Y(x) defined on
[x0-a,x0+a] satisfying Y’(x)=f(x,Y(x)) for xÎ[x0-a,x0+a] and Y(x0)=Y0.
Stability
Y(x0)= Y0+e Y(x) should not vary drastically
Direction Fields
Numerical Solutions
Euler’s Method Y ‘(x)=f(x,Y(x)) xÎ[x0,b] and Y(x0)= Y0
xn=x0+nh h=(xN-x0)/N Note that xN=b
Recall that Y ‘(x)~[Y(x+h)-Y(x)]/h à [Y(xn+1)-Y(xn)]/h ~ f(xn,Y(xn))
à yn+1=yn+hf(xn,yn) n=0,1,…,N-1 y0=Y0
Note that tangent line has slope f(xn,yn) at xn
Example: Y ‘(x)=-Y(x) and Y(0)=1 à true Y(x)=e-x
Euler à yn+1=yn-hyn n=0,…
h=0.1 à y1=y0-hy0=1-(0.1)(1)=0.9 for x1=0.1 error = 0.004837
y2=y1-hy1=0.9-(0.1)(0.9)=0.81 for x2=0.2
error=0.001873
Example: Y ‘(x)=[Y(x)+x2-2]/(x+1) and Y(0)=2
à true Y(x)=x2+2x+2-2(x+1)log(x+1)
Euler à yn+1=yn+h(yn+xn2-2)/(xn+1) n=0,…
Convergence Stability Implicit Methods
Backward
Euler
Y ‘(x)~[Y(x)-Y(x-h)]/h à yn+1=yn+hf(xn+1,yn+1) n=0,… y0=Y0
Example: Y ‘(x) = lY for x>0 and Y(0)=1 à true Y(x)=elx
Euler yn+1= yn+hl yn=(1+hl)yn n=0,…
Backward Euler yn+1=(1-hl)-1yn à =(1-hl)-n
Note that Euler’s method is explicit, i.e., we get yn+1 directly.
However, Backward Euler is implicit in that we must use rootfinding to get yn+1.
Both converge slowly
Trapezoidal
Method
Y ‘(x)=f(x,Y(x)) à Y(xn+1)~Y(xn)+òxnxn+1f(x,Y(x))dx
~Y(xn)+(h/2)[f(xn,Y(xn))+f(xn+1,Y(xn+1))]
à = yn+1=yn+(h/2)[f(xn,yn)+f(xn+1,yn+1)] n=0,… y0=Y0
implicit à requires estimates
Heun’s method - yn+1=yn+(h/2)[f(xn,yn)+f(xn+1,yn)]
Adams-Bashford Method - yn+1=yn+(h/2)[3f(xn,yn)-f(xn-1,yn-1)]
Taylor
Y(xn+1)~Y(xn)+hY ‘(xn)+(h2/2)Y ‘’(xn) with
truncation error term Tn+1=(h3/6)y ‘’’(xn) for xne[xn,xn+1]
Example: Y ’(x) = -Y(x) + 2cos(x) Y(0)=1 true Y(x)=sin(x)+cos(x)
Y ‘’(x) = -Y ‘(x)-2sin(x) = Y(x) - 2cos(x) - 2sin(x)
(derivative of Y ‘(x) )
Y(xn+1)~Y(xn)+h[-Y(xn)+2cos(xn)]+(h2/2)[Y(xn)-2cos(xn)-2sin(xn)]
à y(xn+1)=y(xn)+h[-yn+2cos(xn)]+(h2/2)[yn-2cos(xn)-2sin(xn)] n³0
Runge-Kutta
yn+1=yn+hF(xn,yn;h) y0=Y0 Think of F as average slope
Order 2
F(x,y;h)= υ1f(x,y)+υ2f(x+ah,y+βhf(x,y))
Need to determine constants
Use Taylor expansion for truncation error and solve
υ2¹0, υ1=1-υ2, a=β=1/(2υ2)
favorite choices are υ2= ½ or ¾ or 1.
Order 4 (RK4)
yk+1=yk + h*[f1+2f2+2f3+f4]/6 where f1=f(xkyk),
f2=f(xk+h/2,yk+(h/2)f1), f3=f(xk+h/2,yk+(h/2)f2), and
f4=f(xk+h,yk+hf3)
(Use Taylor expansion as before for order 4)
RKF45 – Runge-Kutta-Fehlberg
yk+1=yk+(16/135)k1+(6656/12825)k3+(28561/56430)k4-(9/50)k5+(2/55)k6
where = k1=hf(xkyk), k2=hf(xk+h/4,yk+k1/4),
k3=hf(xk+3h/8,yk+3k1/32+9k2/32), k4=hf(xk+12h/13,yk+1932k1/2197-7200k2/2197+7296 k3/2197)
k5=hf(xk+h,yk+439k1/216-8k2+3680k3/513-845/4104k4)
k6=hf(xk+h/2,yk-8k1/27+2 2-3544k3/2565+1859k4/4104-11k5/40)
Matlab
ode45 AbsTol RelTol`