Gaussian Elimination

Gaussian Elimination Simultaneous equations, splines, least squares fitting

Simultaneous equations • Gaussian elimination • a11x1+ a12x2 + a13x3 + … a1nxn = b1 • a21x1+ a22x2 + a23x3 + … a2nxn = b2 • a31x1+ a32x2 + a33x3 + … a3nxn = b3 • an1x1 + an2x2 + an3x3 + … annxn = bn • With a set of simultaneous equations, the results are not changed under the operations of • multiplying a row by a scalar • replacing a row by the sum or difference of two rows • interchanging two rows

Pivot elements • Multiply each row i by a11 / ai1 and subtract row 1 from it. • a11x1 + a12x2 + a13x3 + … a1nxn = b1 • a’22x2 + a’23x3 + … a’2nxn = b’2 • etc. • The a11, a’22, a’’33 are the pivot elements • Things to watch for: • if a pivot element is zero, rotate another equation into its place • if a pivot is very small, rotate another equation into its place • in general, use the largest aii as pivot to avoid rounding error. • Numerical Recipes: the Art of Scientific Computing. W. H. Press et al • LINPACK

Electrical Circuit 11 v1 - 5v2 - v6 = 5*V -20 v1 + 41v2 -15v3 - 6v5 = 0 -3v2 + 7v3 - 4v4 = 0 -v3 + 2v4 - v5 = 0 -3v2 - 10v4 + 28v5 - 15v6 = 0 -2v1 - 15v5 + 47v6 = 0

Matrix input MM = 11 -5 0 0 0 -1 -20 41 -15 0 -6 0 0 -3 7 -4 0 0 0 0 -1 2 -1 0 0 -3 0 -10 28 -15 -2 0 0 0 -15 47 >> VV = [50 0 0 0 0 0]

splines

Data and segments Data points: x0,y0 x1,y1 x2,y2xn,yn So, n+1 points The segments are P1 P2Pn Data points: x0,y0 x1,y1 x2,y2 xn-1,yn-1xn,yn The segments will be expressed in terms of second derivatives at the beginning and end of each segment: The segments are P1 P2Pn 2ndderivatives: y”0 y”1 y”2 y”n-1y”n But for a natural spline, the end 2ndderivatives are equal to 0, so there are n-1 to solve for, namely y”1 y”2 y”n-1

Equations to solve 2(x2-x0) (x2-x1) 0 0 0 … 0 0 (x2-x1) 2(x3-x1) (x3-x2) 0 0 … 0 0 0 (x3-x2) 2(x4 – x2) (x4-x3) 0 … 0 0 … 0 0 0 0 0 … (xn-1-xn-2) 2(xn – xn-2) The Y matrix is the column vector of unknowns y”1 y”2 y”3 … y”n-2 y”n-1

example TIME = 0 20 40 60 80 100 120 140 160 180 200 Penicillin = [0 106 1600 3000 5810 8600 9430 10950 10280 9620 9400]

Raw plot

Spline curve >> yy = spline(TIME,Penicillin,XX); >> plot(XX,yy),grid http://www.math.ucla.edu/~baker/java/hoefer/Spline.htm

Least squares fit Zi = A0 + A1 Xi + A2 Xi2 + A3Xi3 + + AnXin What we have is m data points Xi, Yi and we want to choose the A0, A1, etc such that Zi is the best fit to the Yi. What one typically does is minimize the sum of the squares of the distance of Zi from Yi. minimize E = ∑{ Zi - Yi }2 = ∑m { A0 + A1 Xi + A2 Xi2 + A3Xi3 + + AnXin - Yi }2 So, we need partial derivatives wrt each Ai to equal 0.

A0 wrt A0 take derivative of f2 wrtf,then derivative of f wrt A0 = 1 0 = 2 ∑m { A0 + A1 Xi + A2 Xi2 + A3Xi3 + + AnXin - Yi } or ∑m { A0 + A1 Xi + A2 Xi2 + A3 Xi3 + + AnXin } = ∑m Yi or mA0 + A1 ∑m Xi + A2 ∑m Xi2 + A3 ∑m Xi3 + + An ∑mXin = ∑m Yi To do a similar computation wrtA1, take derivative of f2, then derivative of f wrt A1 = Xi

Equations to solve So, finally, the set of simultaneous equations one is solving is: mA0 +A1∑X +A2∑X2 +... +An∑Xn = ∑Y A0∑X +A1∑X2 +A2∑X3 +... +An∑Xn+1 = ∑XY A0∑X2 +A1∑X3 +A2∑X4 +... +An∑Xn+2 = ∑X2Y … A0∑Xn +A1∑Xn+1 + ... +An∑X2n = ∑XnY

example x = -3:.2:2 y = x.^3 + x - 3; plot(x,y),grid

Add error to ‘data’ yy = y + 3*(rand(1,length(x)) - .5); plot(x,y,'-',x,yy,'o'),grid

Best fit [p s] = polyfit(x,yy,3) p = 0.9737 -0.0009 1.1948 -2.5593 yyy = polyval(p,x); % x could be a different vector than the data values plot(x,y,'-',x,yy,'o',x,yyy,'r'),grid

a polyfit version function [ a ] = myPoly(x,y,n) %n contains polynomial order x = (x(:).')'; y = (y(:).')'; k = length(x); xp = 0:1:2*n; xmat = zeros(k,2*n); xmat = [ ones(k,1) xmat]; ymat = repmat(y,1,n+1); for col = 1:2*n xmat(:,col+1) = xmat(:,col).*x; end for col = 2:n+1 ymat(:,col) = ymat(:,col-1).*x; end xvec = sum(xmat); yvec = sum(ymat); MM = zeros(n+1,n+1); C = yvec'; for row = 1:n+1 MM(row,:) = xvec(row:row+n); end a = MM\C; end

Gaussian Elimination