11.27

As described in Sec. PT3.1.2, linear algebraic equations can arise in the solution of differential equations. For example, the following differential equation results from a steady state mass balance for a chemical in a one-dimensional canal.

0 = D * d^c/dx^2 - U * dc/dx - kc

 where c = concentration, t = time, x = distance, D = diffusion coefficient, U = fluid velocity, and k = a first-order decay rate.Convert this differential equation to an equivalent system of simultaneous algebraic equations. Given D = 2, U = 1, k = .2, c(0) = 80 and c(10) = 20, solve these equations from x = 0 to 10 with delta_x = 2, and develop a plot of concentration versus distance.

derivatives

dc_dx = (c(x + dx) - c(x)) ./ dx ;

d2c_dx2 = (c(x + dx) - 2 * c(x) +c(x - dx)) ./ dx^2;

delta x = 1

dc_dx = (c(x + dx) - c(x)) ;

d2c_dx2 = (c(x + dx) - 2 * c(x) +c(x - dx));

0 = D * c(x + dx) - 2 * c(x) + c(x - dx) - U * (c(x + dx) - c(x)

c1 0 = D*c(2) - 2*D*c(1)+ D*c(0) - U*c(2) + U * c(1)

c2 0 = D*c(3) - 2*D*c(2)+ D*c(1) - U*c(3) + U * c(2)

c3 0 = D*c(4) - 2*D*c(3)+ D*c(2) - U*c(4) + U * c(3)

c4 0 = D*c(5) - 2*D*c(4)+ D*c(3) - U*c(5) + U * c(4) 

c5 0 = D*c(6) - 2*D*c(5)+ D*c(4) - U*c(6) + U * c(5)

c6 0 = D*c(7) - 2*D*c(6)+ D*c(5) - U*c(7) + U * c(6)

c7 0 = D*c(8) - 2*D*c(7)+ D*c(6) - U*c(8) + U * c(7)

c8 0 = D*c(9) - 2*D*c(8)+ D*c(7) - U*c(9) + U * c(8)

c9 0 = D*c(10) - 2*D*c(9)+ D*c(8) - U*c(10) + U * c(9)

c1 0 = D*c(2) - U*c(2) + U * c(1) - 2*D*c(1) + D*c(0)

c2 0 = D*c(3) - U*c(3) + U * c(2) - 2*D*c(2) + D*c(1)

c3 0 = D*c(4) - U*c(4) + U * c(3) - 2*D*c(3) + D*c(2)

c4 0 = D*c(5) - U*c(5) + U * c(4) - 2*D*c(4) + D*c(3)

c5 0 = D*c(6) - U*c(6) + U * c(5) - 2*D*c(5) + D*c(4)

c6 0 = D*c(7) - U*c(7) + U * c(6) - 2*D*c(6) + D*c(5)

c7 0 = D*c(8) - U*c(8) + U * c(7) - 2*D*c(7) + D*c(6)

c8 0 = D*c(9) - U*c(9) + U * c(8) - 2*D*c(8) + D*c(7)

c9 0 = D*c(10)- U*c(10)+ U * c(9) - 2*D*c(9) + D*c(8)

Matlab code

D = 2;

U = 1;

k = .2;

%c(0) = 80;

%c(10) = 20;

x = 0:10 ;

dx = 1;% using a delta 1

%        1       2      3     4     5    6    7     8      9

A = [ (U-2*D) (D-U)     0     0     0    0    0     0      0;...

         D   (U-2*D)  (D-U)   0     0    0    0     0      0; ...

         0      D    (U-2*D) (D-U)  0    0    0     0      0; ...

         0      0      D   (U-2*D)  (D-U) 0   0     0      0; ...

         0      0      0     D   (U-2*D)  (D-U) 0   0      0; ...

         0      0      0     0      D   (U-2*D)  (D-U) 0   0; ... 

         0      0      0     0      0   D  (U-2*D) (D-U)   0; ... 

         0      0      0     0      0   0   D  (U-2*D) (D-U); ...

         0      0      0     0      0   0   0     D  (U-2*D);];

B = [-D * 80; 0; 0; 0; 0;0 ; 0 ;0 ; -D*20 + U*20;];

C = A\B;

display(A)

display(B)

display(C)

plot(C)

Results

A =

    -3     1     0     0     0     0     0     0     0

     2    -3     1     0     0     0     0     0     0

     0     2    -3     1     0     0     0     0     0

     0     0     2    -3     1     0     0     0     0

     0     0     0     2    -3     1     0     0     0

     0     0     0     0     2    -3     1     0     0

     0     0     0     0     0     2    -3     1     0

     0     0     0     0     0     0     2    -3     1

     0     0     0     0     0     0     0     2    -3

B =

  -160

     0

     0

     0

     0

     0

     0

     0

   -20

Concentration values =

   79.9413

   79.8240

   79.5894

   79.1202

   78.1818

   76.3050

   72.5513

   65.0440

   50.0293

11.22

Write the following set of equations in matrix form:

50 = 5x_3 -7x_2
4x_2 + 7x_3 + 30 = 0
x_1 - 7x_3 = 40 - 3x_2 +5x_1

 Use excel, MATLAB, or Mathcad to solve for the unkowns. In addition,
 compute the transpose and the inverse of the coefficient matrix.

50 = 5x_3 -7x_2 + 0x_1
30 = - 7x_3 - 4x_2 + 0x_1
-40 =  7x_3 - 3x_2 +4x_1




Matlab Code:

A = [5 -7 0; -7 -4 0; 7 -3 4];

B = [50; 30; -40];

X = A\B;

%transpose

TransposeA = A.';

%inverse

inverseA = inv(A);

display(A)

display(B)

display(X)

display(TransposeA)

display(inverseA)

Results

A =

     5    -7     0

    -7    -4     0

     7    -3     4

B =

    50

    30

   -40

X =     

   -0.1449

   -7.2464

  -15.1812

TransposeA =

     5    -7     7

    -7    -4    -3

     0     0     4

inverseA =

    0.0580   -0.1014         0

   -0.1014   -0.0725         0

   -0.1775    0.1232    0.2500

X_1 = -0.1449

X_2 = -7.2464

X_3 = -15.1812

11.15

Of the following three sets of linear equations, identify the set(s) that you could not solve using an iterative method such as Gauss-seidel. Show using any number of iterations that is necessary that your solution does not converge. Clearly state your convergence criteria (how you know it is not converging).

set one

9x + 3y +z = 13

-6x + 8z =2

2x +5y -z = 6

set two

x +y +6z = 8

x+ 5y - z = 5

4x + 2y - 2z = 4

set three

-3x +4y +5z = 6

-2x +2y -4z = -3

2y -z = 1

Matlab code:

A = [9 3 1; -6 0 8; 2 5 -1];

B = [13; 2; 6];

solution1 = A \ B;

x = [0 0 0];

for ii = 1:100

    x(1) = (B(1) - A(1,2) * x(2) - A(1,3) * x(3))/A(1,1);

    x(2) = (B(2) - A(2,1) * x(1) - A(2,3) * x(3))/A(2,2);

    x(3) = (B(3) - A(3,1) * x(1) - A(3,2) * x(2))/A(3,3);

end

display(solution1);

display(x);

%set two

A = [1 1 6; 1 5 -1; 4 2 -2];

B = [8;5;4];

solution2 =A\B;

x = [0 0 0];

for ii = 1:100

    x(1) = (B(1) - A(1,2) * x(2) - A(1,3) * x(3))/A(1,1);

    x(2) = (B(2) - A(2,1) * x(1) - A(2,3) * x(3))/A(2,2);

    x(3) = (B(3) - A(3,1) * x(1) - A(3,2) * x(2))/A(3,3);

end

display(solution2);

display(x);

%set three

A = [-3 4 5; -2 2 -4; 0 2 -1];

B = [6; -3; 1];

solution3 =A\B;

x = [0 0 0];

for ii = 1:100

    x(1) = (B(1) - A(1,2) * x(2) - A(1,3) * x(3))/A(1,1);

    x(2) = (B(2) - A(2,1) * x(1) - A(2,3) * x(3))/A(2,2);

    x(3) = (B(3) - A(3,1) * x(1) - A(3,2) * x(2))/A(3,3);

end

display(solution3);

display(x);

Results

solution1 =

   1.0000e+00

   1.0000e+00

   1.0000e+00

x =

  -Inf  -Inf  -Inf

solution2 =

     1

     1

     1

x =

 -2.3612e+101  5.5271e+100 -4.1696e+101

solution3 =

   6.9565e-01

   9.3478e-01

   8.6957e-01

x =

  -1.6803e+93  -3.1205e+93  -6.2411e+93

All 3 sets diverge

10.25

Polynomial interpolation consists of determining the unique (n - 1)th-order polynomial that fits n data points. Such polynomials have the general form.

 f(x) = P1x^n-1 + P2x^n-2 +... + Pn-1x +Pn                (eq10.25)

 where the p's are constant coefficients. A straightforward way for computing the coefficients is to generate n linear algebraic equations that we can solve simultaneously for the coefficients. Suppose that we want to determine the coefficients of the fourth-order polynomial f(x) =  p1x^4 + p2x^3 +p3x^2 +p4x +p5 that passes through the following five points:(200,0.746),(250, .675), (300,.616), (400, .525) and (500, .457). Each of these pairs can be substituted into Eq. (P10.25) to yield a system of five equations with five unknowns (the p's). Use this approach to solve for the coefficients. In addition, determine and interpret the condition number.


Matlab code:

B = [ .746; .675; .616; .525; .457]

x = [200, 250, 300, 400, 500];

xx = [1,1,1,1,1];

A = [ x.^4 ;x.^3; x.^2 ;x.^1; xx]'

p = A\B

A * p

cond(A)

RESULTS

B =

   7.4600e-01

   6.7500e-01

   6.1600e-01

   5.2500e-01

   4.5700e-01

A =

   1.6000e+09   8.0000e+06   4.0000e+04   2.0000e+02   1.0000e+00

   3.9062e+09   1.5625e+07   6.2500e+04   2.5000e+02   1.0000e+00

   8.1000e+09   2.7000e+07   9.0000e+04   3.0000e+02   1.0000e+00

   2.5600e+10   6.4000e+07   1.6000e+05   4.0000e+02   1.0000e+00

   6.2500e+10   1.2500e+08   2.5000e+05   5.0000e+02   1.0000e+00

p =

   1.3333e-12                P1

  -4.5333e-09                P2

   5.2967e-06                P3

  -3.1737e-03                P4

   1.2030e+00               P5

ans =

   7.4600e-01

   6.7500e-01

   6.1600e-01

   5.2500e-01

   4.5700e-01

cond(A) =

   1.1712e+13

A has a large condition number, meaning it is ill-formed, and there is a high potential for error

10.21

Consider vectors:
 A = 2i - 3j + ak
 B = bi + j - 4k
 C = 3i + cj +2k

 Vector A is perpendicular to B as well as to C. It is also known that B * C = 2. Use any method studied in this chapter to solve for the three unknowns, a, b, and c.


Matlab code:

% A * B = 0 ... 2bi -3j - 4ak = 0 ... 3j

% A * C = 0 ... 6i - 3cj + 2ak = 0... 6i

% B * C = 2 ... 3bi + cj - 8k = 2k...  10k

%       a      b      c

A = [-4 2 0; 2 0 -3; 0 3 1]

b = [3; -6; 10]

x = A\b

display(A*x)

First take the dot products. Since A is perpendicular to B as well as to C we can say that they are orthogonal so their dot products are 0.  We then pull the vectors 3j, 6i, and 8k from the dot products, leaving us with a b matrix of [3;6;10]. then solves for the unknowns, giving us A=0.5250, B=2.5500, and C=2.3500.  We confirm the result by multiplying A and x, reproducing b

A =

    -4     2     0

     2     0    -3

     0     3     1

b =

     3

    -6

    10

x =

    0.5250

    2.5500

    2.3500

ans =

    3.0000

   -6.0000

   10.0000

10.15

(a)Determine the condition number for the following system using
 the row-sum norm. Do not normalize the system.

A = [ 1 4 9 16 25;4 9 16 25 36; 9 16 25 36 49; 16 25 36 49 64;25 36 49 64 81]

How many digits of precision will be lost due to ill-conditioning?
(b)Repeat (a), but scale the matrix by making the maximum element in each
row equal to one.


Matlab Code:

A = [ 1 4 9 16 25;4 9 16 25 36; 9 16 25 36 49; 16 25 36 49 64;25 36 49 64 81]

%(a)

conditionA = cond(A)

significant = 10^(log10(cond(A)) - 7.2)

%(b)

b = [ 1/25 0 0 0 0;0 1/36 0 0 0;0 0 1/49 0 0;0 0 0 1/64 0;0 0 0 0 1/81];

A_scaled = b * A

conditionA_scaled = cond(A_scaled)

significant_scaled = 10^(log10(cond(A_scaled)) - 7.2)

RESULTS

A =

     1     4     9    16    25

     4     9    16    25    36

     9    16    25    36    49

    16    25    36    49    64

    25    36    49    64    81

conditionA =

     3.145210075463109e+17

significant =

     1.984493397046558e+10

A_scaled =

   0.040000000000000   0.160000000000000   0.360000000000000   0.640000000000000   1.000000000000000

   0.111111111111111   0.250000000000000   0.444444444444444   0.694444444444444   1.000000000000000

   0.183673469387755   0.326530612244898   0.510204081632653   0.734693877551020   1.000000000000000

   0.250000000000000   0.390625000000000   0.562500000000000   0.765625000000000   1.000000000000000

   0.308641975308642   0.444444444444444   0.604938271604938   0.790123456790123   1.000000000000000

conditionA_scaled =

     9.237149093911445e+16

significant_scaled =

     5.828247062861999e+09

These results show that both matrices are substantially ill-formed and there is a large loss of significance

10.8

The following system of equations is designed to determine concentrations (the c's in g/m^3) in a series of coupled reactors as a function of the amount of mass input to each reactor (the right-hand sides in g/day),

15c1 - 3c2 - c3 = 3300
-3c1 +18c2 - 6c3 = 1200
-4c1 - c2 + 12c3 = 2400

(a) Determine the matrix inverse(b) Use the inverse matrix to determine a solution(c) Determine how much the rate of mass input to reactor 3 must be increased to include a 10 g/m^3 rise in the concentration of reactor 1.(d)How much will the concentration in reactor 3 be reduced if the rate of mass input to reactors 1 and 2 is reduced by 700 and 350 g/day, respectively?


Matlab code:

%(a) Determine the matrix inverse

A = [15 -3 -1; -3 18 -6; -4 -1 12]

inversematrix = inv(A)

%(b) Use the inverse matrix to determine a solution

B = [3300; 1200; 2400]

x = inversematrix * B

OriginalB = A * x

%(c) Determine how much the rate of mass input to reactor 

%3 must be increased to incduce a 10 g/m^3 rise in the concentration of

%reactor 1.

[x(1)+10; x(2); x(3)]

OriginalBplus10 = A * [x(1)+10; x(2); x(3)]

%(d)How much will the concentration in reactor 3 be reduced if the rate of

%mass input to reactors 1 and 2 is reduced by 700 and 350 g/day,

%respectively?

Bnew = [2600; 850; 2400]

newx = inversematrix * Bnew

newB = A * newx

diffx = x - newx

Results

(A)

A =

    15    -3    -1

    -3    18    -6

    -4    -1    12

inversematrix =

   0.072538860103627   0.012780656303972   0.012435233160622

   0.020725388601036   0.060794473229706   0.032124352331606

   0.025906735751295   0.009326424870466   0.090155440414508

(B)

B =

        3300

        1200

        2400

x =

   1.0e+02 *

   2.845595854922280

   2.184455958549223

   3.130569948186529

OriginalB =

   1.0e+03 *

   3.299999999999999

   1.200000000000001

   2.400000000000000

ans =

   1.0e+02 *

   2.945595854922280

   2.184455958549223

   3.130569948186529

(C)

OriginalBplus10 =

   1.0e+03 *

   3.450000000000000

   1.170000000000001

   2.360000000000000

(D)

Bnew =

        2600

         850

        2400

newx =

   1.0e+02 *

   2.293091537132988

   1.826597582037997

   2.916580310880829

newB =

   1.0e+03 *

   2.600000000000000

   0.850000000000001

   2.400000000000000

diffx =

  55.250431778929169

  35.785837651122620

  21.398963730569960

10.5

Determine the total flops as a function of the number of equations n for the  (a) decomposition, (b) forward-substitution, and (c) back-substitution phases of the LU decomposition version of Gauss elimination.

Decomposition
n^3/3 - n/3 or as n increases n^3/3 + O(n)

Forward-substitution
n^3/3 - n/3 + n^2 or as n increases n^3/3 + O(n^2) 

Back-substitution
n^3/3 - n/3 + n^2 or as n increases n^3/3 + O(n^2)