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)e^j2Pftdf and e^jq=cos(q)+j*sin(q)

Discrete: x₀, x₁,…, x_N-1 with x_i – x_real+j* x_imag

X(n) = (1/n)*S_k=1^N-1 x(k)*e^-jk2Pn/N Note that x(n) = S_k=1^N-1 X(k)*e^jk2Pn/N for n=0,…,N-1

Magnitude = ||X(n)|| = (x_real*x_real + x_imag*x_imag)^0.5

Phase = tan^-1(x_imag/x_real)

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 – X_k à (1/N)x_N-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) = lim_hà0 [f(x+h)-f(x)]/h

f ’(x) ~ [f(x+h)-f(x)]/h = D_hf(x) numerical derivative for small h (stepsize)

forward difference

Taylor series à f(x+h) = f(x) + hf ’(x) + (h²/2)*f ‘’(c) for cÎ[x,x+h]

à D_hf(x) = (1/h)*[f(x)+hf ’(x)+(h²/2)*f’ ’’(c)-f(x)] = f ’(x)+(h/2)*f ‘’(c)

à error = f ’(x) - D_hf(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

D_hf(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, x₀ = x₁-h, x₁, x₂ = x₁+h y_i=f(x_i) for i=0,1,2

P₂(x)=y₀*(x-x₁)(x-x₂)/(x₀-x₁)(x₀-x₂)+y₁*(x-x₀)(x-x₂)/(x₁-x₀)(x₁-x₂)+

y₂*(x-x₀)(x-x₁)/(x₂-x₀)(x₂-x₁)

= y₀*(x-x₁)(x-x₂)/2h²+y₁*(x-x₀)(x-x₂)/(-h²)+y₀*(x-x₀)(x-x₁)/(2h²)

P₂’(x) = (2x-x₁-x₂)y₀/2h²+(2x-x₀-x₂)y₁/(-h²)+(2x-x₀-x₁)y₂/2h²

à P₂’(x₁) = (x₁-x₂)y₀/2h²+(2x₁-x₀-x₂)y₁/(-h²)+(x₁-x₀)y₂/2h² = [f(x₂)-f(x₀)]/2h

à f’(x₁) ~ P₂’(x) = [f(x₁+h)-f(x₁-h)]/2h = D_hf(x₁) Central difference formula

f’(x)-P_n’(x) = Y_n(x)*f⁽ⁿ⁺²⁾(c₁)/(n+2)! + Y_n’(x)*f⁽ⁿ⁺¹⁾(c₂)/(n+2)! With Y_n(x)=P_j=1ⁿ(x-x_j)

c₁ and c₂Î[min x_j value, max x_j value] and the x_j values include x for this interval

n=2 yields Y_n(x)=(x-x₀)(x-x₁)(x-x₂) and

Y_n’(x)=(x-x₁)(x-x₂)+(x-x₀)(x-x₂)+(x-x₀) (x-x₁)

à Y_n’(x₁)=(x₁-x₀)(x₁-x₂) = -h²

à f’(x₁) ~ [f(x₁+h)-f(x₁-h)]/2h = (h²/6)*f”’(c₂)

Undetermined Coefficients

f ‘’(x)~D_h⁽²⁾f(x)=Af(x+h)+Bf(x)+Cf(x-h)

Taylor à f(x-h)~f(x)-hf ’(x)+h²/2*f ‘’(x)-h³/6*f ‘’’(x)+h⁴/24*f ‘’’’(x)

à f(x+h)~f(x)+hf ’(x)+h²/2*f ‘’(x)+h³/6*f ‘’’(x)+h⁴/24*f ‘’’’(x)

à D_h⁽²⁾f(x) ~ (A+B+C)*f(x) + h*(A-C)*f ’(x) h²/2*(A+C)*f ‘’(x) +

h³/6*(A-C)*f ‘’’(x)+h⁴/24*(A+C)*f ‘’’’(x)

D_h⁽²⁾f(x) ~ f ”(x) à (A+B+C) = coefficient of f(x) = 0;

h(A-C) = coefficient of f ’(x) = 0;

h²/2*(A+C) = coefficient of f ‘’(x) = 1

à A = C = 1/h² and B = -2/h²

à D_h⁽²⁾f(x) ~ [f(x+h)-2f(x)+f(x-h)]/h²

à f ‘’(x)- D_h⁽²⁾f(x) ~ -h²/12*f ‘’’’(x)

Numerical Integration (Quadrature)

I(f) = ò_a^bf(x)dx = F(b)-F(a) where F(x) is antiderivative of f(x)

Newton-Cotes Formulae

I(f) = ò_a^bf(x)dx ~ ò_a^bP_n(x)dx where P_n(x) = å_j=0ⁿa_jx^j = a polynomial

Trapezoidal Rule

P₁(x) = [(b-x)f(a)+(x-a)f(b)]/(b-a) -->

T₁(f) = ò_a^b[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)³

Example: f(x) = 1/(1+x) a=0 and b=1 Note: F(x) = ln(1+x)

T₁ = 0.75 but I(f) = ln(2) = 0.693147

Composite: n intervals h=(b-a)/n x_j=a+h*j j=0,1,…,n

I(f) = ò_aj^b f(x)dx =å_j=0^n-1ò_xj^xj+1 f(x)dx ~ Tn(f) = (h/2)å_j=0^n-1[f(x_j)+f(x_j+1)]

= (h/2)[f(x₀)+f(x_n)+2å_j=1^n-1f(xj)]=(b-a)*[f(x₀)+f(x_n)+2å_j=1^n-1f(xj)]/2n

where (b-a) = width and [f(x₀)+f(x_n)+2å_j=1^n-1f(xj)]/2n = average heigth

Error -(b-a)³å_j=1ⁿ f ‘’(c_j)/12n³

Example: n = 4

h=(b-a)/4 x₀=a, x₁=a+h=(3a+b)/4, x₂=a+2h=(a+b)/2, x₃=a+3h=(a+3b)/4, x₄=a+4h=b

f(x) = e^(-x²), n=4, a=0, b=1 --> T₄=0.7468241…

f(x) = 1/(1+x²), n=4, a=0, b=4 --> T₄=1.32581766…

I = tan^-1(4)

f(x) = 1/(2+cos(x)), n=4, a=0, b=2p --> T₄=0.7468241…

I=2p/Ö3

Simpson’s Rule

N=2 c=(a+b)/2 h=(b-a)/n

I(f) ~ ò_a^bP₂(x)dx = ò_a^b [(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

=ò_a^b [(x-c)(x-b)f)a)/(a-c)(a-b)]dx+ò_a^b[(x-a)(x-b)f(c)/(a-c)(b-c)]dx+

ò_a^b [(x-a)(x-c)f(b)/(b-a)(b-c)]dx

Note: ò_a^b(x-c)(x-b)/(a-c)(a-b)dx = (1/2h²)ò_a^a+nh(x-c)(x-b)dx

=(1/2h²)ò_a^2h(u-h)(u-2h)du with u=x-a

=(1/2h²)[u³/3-3u²h/2+2h²u]₀^2h = h/3

à I(f) ~ S₂(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

S₂(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 x_j=a+h*j j=0,1,…,n

I(f)=ò_a^bf(x)dx=ò_x0^xn f(x)dx=ò_x0^x2f(x)dx +ò_x2^x4f(x)dx +…+ò_xn-2^xn f(x)dx

S_n(f) = (h/3)[f(x₀)+4f(x₁)+2f(x₂)+4f(x₃)+2f(x₄)+…+2f(x_n-2)+4f(x_n-1)+f(x_n)]

=(b-a)[f(x₀)+f(x_n)+4å_{j=1,j odd}^n-1f(x_j)+2å_{j=1,j even}^n-1f(xj)]/6

Error = [(b-a)⁵ /180n⁴]f⁽⁴⁾ where f⁽⁴ = average f⁽⁴⁾ value

Composite’ use of 3^rd order Lagrange polynomial

I(f) ~ (3h/8) [f(x₀)+3f(x₁)+3f(x₂)+ f(x₃]

=(b-a) [f(x₀)+3f(x₁)+3f(x₂)+ f(x₃]/8

Error = -[3h⁵/80]*f⁽⁴⁾(c) = -[(b-a)⁵/6480]*f⁽⁴⁾(c)

Can do higher orders

Gaussian Quadrature

Assume a = -1 and b = 1 I(f) = ò_a^bf(x)dx = ò_-1¹f(x)dx ~ I_n(f) = S_j=1nw_jf(x_j)

Choose x and w values so I_n(f) = I(f) for as large a degree n as possible for simple functions

n = 1

f(x) =1 à I(f) = ò_-1¹dx = 2 = I₁(f) = w₁ à w₁ = 2

f(x) = x à I(f) = ò_-1¹xdx = 0 = I₁(f) = w₁x₁ = 2x₁ à x₁ = 0

à I₁(f) = 2f(0)

Note: I(f) = ò_a^bf(x)dx ~ ò_a^bf[(a+b)/2]dx = (b-a)f[(b+a)/2]

a = -1 and b = 1 à I(f) ~ 2f(0) midpoint formula

n = 2

I₂(f) = w₁f(x₁)+w₂f(x₂)

f(x)=1,x,x²,x³ à w₁+w₂=2; w₁x₁+w₂x₂=w₁x₁³+w₂x₂³=0; w₁x₁²+w₂x₂²=2/3

à w₁=w₂=1, x₁=-Ö3/3, x₂=Ö3/3 à I₂(f) = f(-Ö3/3)+f(Ö3/3)

n>2

I_n(f) = S_j=1ⁿw_jf(x_j)

f(x)= 1,x,x²,…,x^2n-1 à S_j=1ⁿ w_j=0=S_j=1ⁿ w_jfx_j³, 2/3=S_j=1ⁿ w_jx_j², … ,

S_j=1ⁿ w_jx_j^2n-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ò_-1¹f([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) = ò₀¹(1/Öx)dx = 2 = I₁(f) = w₁ à w₁ = 2

f(x) = x à I(f) = ò₀¹(x/Öx)dx = 2/3 = I₁(f) = w₁x₁ = 2x₁ à x₁ = 1/3

à I₁(f) = 2f(1/3)

n = 2

I₂(f) = w₁f(x₁)+w₂f(x₂)

f(x)=1,x,x²,x³ à w₁+w₂=2; w₁x₁+w₂x₂=2/3; w₁x₁²+w₂x₂²=2/5; w₁x₁³+w₂x₂³=2/7

à w₁=1+Ö30/18; w₂=1-Ö30/18; x₁=3/7-2Ö30/35; x₂=3/7+2Ö30/25

n>2

S_j=1ⁿ w_j=2; S_j=1ⁿ w_jfx_j=2/3; S_j=1ⁿ w_jx_j²=2/5; … ; S_j=1ⁿ w_jx_j^2n-1=2/(4n-1)

Linear Equations

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=1ⁿa_ijx_j=b_i 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 a_ij = max(i,j), b_i=1 for all i,j à x_j =0 for j=1,..,n-1 and x_n=1/n. If n = 3, we have a₁₁=1, a₁₂=2, a₁₃=3, a₁₂=a₂₂=2, a₂₃=3, a₁₃=a₂₃=a₃₃=3, b₁=b₂=b₃=1 à x₁=x₂=0 and x₃=1/3.

Linear Algebra - Matrix arithmetic

Scalar – number

Vector – one-dimensional array of numbers

Matrix – two-dimensional array of numbers

Matrix A_mXn

Transpose B_nXm=A_nXm^T à b_ij= a_ji for all i,j

Scalar multiplication - B_mXn=sA_mXn à b_ij= s*a_ij for all i and j

Addition (Subtraction) - C_mXn=A_mXn±B_mXn à c_ij=a_ij±b_ij for all i and j

Multiplication - C_mXs=A_mXn±B_rXs à c_ij=å_k=1ⁿa_ik*b_kj for all i and j and where n = r

Note: c_ij = dot or vector product of i^th row of A (a vector) and j^th column of B (a vector)

Back to equations:

Let A = matrix with m rows and n columns = {a_ij}_mCn, b = (column) vector with m rows = {b_i}_mC1, and x = (column) vector with n rows = {x_j}_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 {a_ii} 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, I_ij=δ_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=B^TA^T; (A+B)^T=A^T+B^T; (cA)^-1=(1/c)A^-1; (AB)^-1=B^-1A^-1; det(AB)=det(A)det(B); det(A^T)=det(A); det(cA)=cⁿdet(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) = a₁₁a₂₂- a₂₁a₁₂a=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

More secrets

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 x₁ 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 x₂ 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 ½* x₃=1/2, then x₃=1. So, equation 2 gives us –2x₂+ x₃=-2x₂+1=3 or x₂=-1. This, in turn, can be used in equation 1 to give x₁+2x₂+x₃= x₁+2(-1)+1= x₁-1=0 or x₁=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=1ⁿa_ij⁽¹⁾x_j=b_i⁽¹⁾ for i=1,…n. Now, for k=1,2,…,n-1 do the elimination steps: i) define multipliers m_ik=a_ik^(k)/a_kk^(k) assuming a_kk^(k)¹0 for i=k+1,…,n; ii) subtract m_ik*E(k) from E(i) to eliminate xk from E(i), yielding a_ij^(k+1)= a_ij^(k)-m_ika_kj^(k) for i,j=k+1,…n and b_i^(k+1)=b_i^(k)-m_ik b_k^(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 u_ij=a_ij⁽ⁿ⁾ and g_i=b_i⁽ⁿ⁾, use back substitution to solve, giving i) x_n= g_n/u_nn; and ii) x_i=[g_i-å_j=i+1ⁿu_ijx_j]/u_ii for i=n-1,…,1.

One problem is what to do if a_kk^(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.

Back to Don’s dirty secret method

Basically, consider the E matrix where E=[A|b|I]. Note that E has m rows and 2n+1 columns. Note that e_ij = a_ij for j=1,…,n; e_in+1 = b_i, and e_ijj=δ_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 e_rs¹0 then we look for e_is¹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 e_rs so e ‘_rj=(1/e_rs)*e_rj 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 -e_is and add it to row i, i.e., e’_ij=(-e_is)*e’_rj+e_ij=(-e_is)*(e_rj/e_rs)+e_ij for j=1,2,…,2n+1. Note that e_is=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=[x₁;x₂], 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 e_rs=e₁₁=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/e_rs. Now, multiplying the new row 1 by the element in the pivot column but in the next row i=2, –e_is=-e₂₁=-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 e_rs=e₂₂=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 –e₁₂=-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 u_ij=a_ij⁽ⁱ⁾ 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 S_j=1ⁿa_ij⁽¹⁾x_j=b_i⁽¹⁾ with m_ij=a_ij^(j)/a_jj^(j) "i=j+1,…,n à a_ij^(k+1)=a_ij^(k)-m_ika_kj^(k) "i,j=k+1,…,n and b_i^(k+1i)=b_i^(k)-m_ikb_k^(k) "i=k+1,…,n

We can define L via l_ij=m_ij⁽ⁱ⁾ 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]

m₂₁=3/4, m₃₁=1/2, m₄₁=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, m₃₂=6/7 and m₄₂=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, m₄₃=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 - a_ij=0 if i¹j

Tridiagonal matrix - a_ij=0 if |i-j|>1

Let a_ii=b_i, a_ii+1=c_i, and a_ii-1=a_i so u_ii=b_i and u_ii+1=c_i and l_ii=1 and l_ii-1=a_i

à b₁=b₁, a_jb_j-1=a_j, a_jc_j-1+b_j=b_j àb₁=b₁, a_j=a_j/b_j-1, b_j=b_j–a_jc_j-1 "j=2,…,n

Consider Ax=f à Lg=f, Ux=g à g₁=f₁, g_j=f_j-a_jg_j-1 "j=2,3,…,n

à x_n=g_n/b_n, x_n-1j=[g_j-c_jx_j+1]/b_j "j=n-1,…,1

Iterative Methods Ax=b

Jacobi Iteration

Ax=b à S_j=1ⁿa_ijx_j=S_j=1,j¹iⁿa_ij⁽¹⁾x_j+a_iix_i=b_i⁽¹⁾b_i à x_i=(1/a_ii)[b_i-S_j=1,j¹1ⁿa_ijx_j] with a_ii¹0

Let x⁰=initial estimate for vector x. Then, x_i^(k+1)=(1/a_ii)[b_i-S_j=1,j¹1ⁿa_ijx_j^(k)]

Gauss-Seidel

Similar to Jacobi but x_i^(k+1)=(1/a_ii)[b_i-S_j=1,j>1ⁿa_ijx_j^(k)--S_j=1,j<1ⁿa_ijx_j^(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.

Determinants

Consider a square matrix A of order n.

det(A) = |A| =å _j=1ⁿa_ijC_ij where C_ij = a cofactor = (-1)^i+jM_ij

where M_ij = 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)

Rules

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=1ⁿa_jj

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={a_ij} and lI-A=[l-a₁₁ -a₁₂; -a₂₁ l-a₂₂] so f(l)=l²-(a₁₁+a₂₂)l+a₁₁a₂₂-a₂₁a₁₂. Let A=[1.25 0.75; 0.75 1.25]. Note det(A)=4. f(l)=l²-2.5l+1=0 à l₁=0.5 and l₂=2 as roots or eigenvalues.

Solving for v^j, the eigenvector (a value of v to solve equation for a given eigenvalue) for eigenvalue l_j, 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] and v⁽²⁾=[-1;1] will work.

Example: n=3. A=[-7 13 -16; 13 -10 13; -16 13 -7]. f(l)=l³+24l²-405l+972=0 so l₁=-36, l₂=9, and l₃=3. Suppose we select l=l₁=-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 v₁⁽¹⁾=1 arbitrarily. Then, -13v₂⁽¹⁾+16v₃⁽¹⁾=29,

-26v₂⁽¹⁾-13v₃⁽¹⁾=13, and -13v₂⁽¹⁾-29v₃⁽¹⁾= -16. Solving, we can get v₂⁽¹⁾=-1 and v₃⁽¹⁾=1.

Symmetric Matrices

For a symmetric matrix, where A^T=A or a_ij=a_ji 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)||=ÖS_j=1ⁿ v^{(r) 2}||=1). Moreover, for any vector x with n elements, there are unique constants c_j such that x=S_j=1ⁿ c_jv^(j) and c_i=S_j=1ⁿ x_jv_j⁽ⁱ⁾ .In addition, if we define U={v^(j)} then D=U^TAU = {l_jδ_ij} and UU^T=U^TU=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 |l₁| > |l₂| ³…³|l_n|. This method calculates l₁ and v⁽¹⁾. Let z⁽⁰⁾=initial guess for v⁽¹⁾, usually done randomly. Define w⁽¹⁾=Az⁽⁰⁾ and a¹=largest component of w⁽¹⁾. Now, iteratively, w^(m)=Az^(m-1), a^m= maximum component of w^(m), z^(m)=w(m)/ a^m. Moreover, l₁^(m)=w_k^(m)/z_k^(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 l₁as mà¥.

Nonlinear equation systems

f(x,y)=0, g(x,y)=0

Newton’s method

(x₀,y₀)=initial guess

r(x,y)=f(x₀,y₀)+(x-x₀)f_x(x₀,y₀)+(y-y₀)f_y(x₀,y₀) where f_x(x,y)=δf(x,y)/δx and f_y(x,y)=δf(x,y)/δy. Note that r(x,y) = equation of tangent plane to z=f(x,y) at (x₀,y₀,f(x₀,y₀)).

q(x,y)=g(x₀,y₀)+(x-x₀)g_x(x₀,y₀)+(y-y₀)g_y(x₀,y₀) where g_x(x,y)=δg(x,y)/δx and g_y(x,y)=δg(x,y)/δy. Note that q(x,y) = equation of tangent plane to z=g(x,y) at (x₀,y₀,g(x₀,y₀)).

Need to solve r(x,y)=0 and q(x,y)=0 à x₁=x₀+δx and y₁=y₀+δy. Repeat iteratively.

Generalized method Modified Method

Ordinary Differential Equations(ODE)

Y’(x)=dy/dx=f(x,Y(x) x³(x₀) 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(x₀)= Y₀

Example: Y ’(x)=sin(x)à Y(x)= - cos(x)+c x₀=P/3 and Y(x₀)=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³ x₀

à d(e^-lxY(x))/dx=e^-lxb(x) à e-xY(x)=c+ò_x0^xe^-ltb(t)dt

à c=e^-lx0Y(x₀)=e^-lx0Y₀

Example: Y ‘(x)=-[Y(x)]²+Y(x) Y₀=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 (x₀,Y₀), then $ unique function Y(x) defined on

[x₀-a,x₀+a] satisfying Y’(x)=f(x,Y(x)) for xÎ[x₀-a,x₀+a] and Y(x₀)=Y₀.

Stability

Y(x₀)= Y₀+e Y(x) should not vary drastically

Direction Fields

Numerical Solutions

Euler’s Method Y ‘(x)=f(x,Y(x)) xÎ[x₀,b] and Y(x₀)= Y₀

x_n=x₀+nh h=(x_N-x₀)/N Note that x_N=b

Recall that Y ‘(x)~[Y(x+h)-Y(x)]/h à [Y(x_n+1)-Y(x_n)]/h ~ f(x_n,Y(x_n))

à y_n+1=y_n+hf(x_n,y_n) n=0,1,…,N-1 y₀=Y₀

Note that tangent line has slope f(x_n,y_n) at x_n

Example: Y ‘(x)=-Y(x) and Y(0)=1 à true Y(x)=e^-x

Euler à y_n+1=y_n-hy_n n=0,…

h=0.1 à y₁=y₀-hy₀=1-(0.1)(1)=0.9 for x₁=0.1 error = 0.004837

y₂=y₁-hy₁=0.9-(0.1)(0.9)=0.81 for x₂=0.2

error=0.001873

Example: Y ‘(x)=[Y(x)+x²-2]/(x+1) and Y(0)=2

à true Y(x)=x²+2x+2-2(x+1)log(x+1)

Euler à y_n+1=y_n+h(y_n+x_n²-2)/(x_n+1) n=0,…

Convergence Stability Implicit Methods

Backward Euler

Y ‘(x)~[Y(x)-Y(x-h)]/h à y_n+1=y_n+hf(x_n+1,y_n+1) n=0,… y₀=Y₀

Example: Y ‘(x) = lY for x>0 and Y(0)=1 à true Y(x)=e^lx

Euler y_n+1= y_n+hl y_n=(1+hl)y_n n=0,…

Backward Euler y_n+1=(1-hl)^-1y_n à =(1-hl)^-n

Note that Euler’s method is explicit, i.e., we get y_n+1 directly.

However, Backward Euler is implicit in that we must use rootfinding to get y_n+1.

Both converge slowly

Trapezoidal Method

Y ‘(x)=f(x,Y(x)) à Y(x_n+1)~Y(x_n)+ò_xn^xn+1f(x,Y(x))dx

~Y(x_n)+(h/2)[f(x_n,Y(x_n))+f(x_n+1,Y(x_n+1))]

à = y_n+1=y_n+(h/2)[f(x_n,y_n)+f(x_n+1,y_n+1)] n=0,… y₀=Y₀

implicit à requires estimates

Heun’s method - y_n+1=y_n+(h/2)[f(x_n,y_n)+f(x_n+1,y_n)]

Adams-Bashford Method - y_n+1=y_n+(h/2)[3f(x_n,y_n)-f(x_n-1,y_n-1)]

Taylor

Y(x_n+1)~Y(x_n)+hY ‘(x_n)+(h²/2)Y ‘’(x_n) with

truncation error term Tn+1=(h3/6)y ‘’’(x_n) for x_ne[x_n,x_n+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(x_n+1)~Y(x_n)+h[-Y(x_n)+2cos(x_n)]+(h²/2)[Y(x_n)-2cos(x_n)-2sin(x_n)]

à y(x_n+1)=y(x_n)+h[-y_n+2cos(x_n)]+(h²/2)[y_n-2cos(x_n)-2sin(x_n)] n³0

Runge-Kutta

y_n+1=y_n+hF(x_n,y_n;h) y₀=Y₀ Think of F as average slope

Order 2

F(x,y;h)= υ₁f(x,y)+υ₂f(x+ah,y+βhf(x,y))

Need to determine constants

Use Taylor expansion for truncation error and solve

υ₂¹0, υ₁=1-υ₂, a=β=1/(2υ₂)

favorite choices are υ₂= ½ or ¾ or 1.

Order 4 (RK4)

y_k+1=y_k + h*[f₁+2f₂+2f₃+f₄]/6 where f₁=f(x_ky_k),

f₂=f(x_k+h/2,y_k+(h/2)f₁), f₃=f(x_k+h/2,y_k+(h/2)f₂), and

f₄=f(x_k+h,y_k+hf₃)

(Use Taylor expansion as before for order 4)

RKF45 – Runge-Kutta-Fehlberg

y_k+1=y_k+(16/135)k₁+(6656/12825)k₃+(28561/56430)k₄-(9/50)k₅+(2/55)k₆

where = k₁=hf(x_ky_k), k₂=hf(x_k+h/4,y_k+k₁/4),

k₃=hf(x_k+3h/8,y_k+3k₁/32+9k₂/32), k₄=hf(x_k+12h/13,y_k+1932k₁/2197-7200k₂/2197+7296 k₃/2197)

k₅=hf(x_k+h,y_k+439k₁/216-8k₂+3680k₃/513-845/4104k₄)

k6=hf(x_k+h/2,y_k-8k₁/27+2 ₂-3544k₃/2565+1859k₄/4104-11k₅/40)

Matlab

ode45 AbsTol RelTol`