1
CHAPTER 9
9.1 The flop counts for the tridiagonal algorithm in Fig. 9.6 can be summarized as
Mult/Div Add/Subtr Total
Forward elimination 3(n – 1) 2(n – 1) 5(n – 1)
Back substitution 2n – 1 n–1 3n – 2
Total 5n – 4 3n – 3 8n – 7
Thus, as n increases, the effort is much, much less than for a full matrix solved with Gauss elimination
which is proportional to n3.
9.2 The equations can be expressed in a format that is compatible with graphing x2 versus x1:
18 3 x1
x2 3 0.5 x1
6
40 2 x1
x2 5 0.25 x1
8
which can be plotted as
a11=3;a12=-6;b1=-18;
a21=-2;a22=8;b2=40;
x1=[0:20];x21=(b1-a11*x1)/a12;x22=(b2-a21*x1)/a22;
plot(x1,x21,x1,x22,'--'),grid
xlabel('x_1'),ylabel('x_2')
14
12
10
x2
2
0 5 10 15 20
x1
Thus, the solution is x1 = 8, x2 = 7. The solution can be checked by substituting it back into the equations to
give
3(8) 6(7) 24 42 18
2(8) 8(7) 16 56 40
9.3 (a) The equations can be rearranged into a format for plotting x2 versus x1:
x2 14.25 0.77 x1
1.2
x2 11.76471 x1
1.7
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
2
20
0
0 20 40 60 80
-20
-40
-60
If you zoom in, it appears that there is a root at about (38.7, –15.6).
38 38.5 39 39.5
-15
-15.4
-15.8
-16.2
The results can be checked by substituting them back into the original equations:
0.77(38.7) 15.6 14.2 14.25
1.2(38.7) 1.7(15.6) 19.92 20
(b) The plot suggests that the system may be ill-conditioned because the slopes are so similar.
(c) The determinant can be computed as
D 0.77(1.7) 1(1.2) 0.11
which is relatively small. Note that if the system is normalized first by dividing each equation by the largest
coefficient,
0.77 x1 x2 14.25
0.705882 x1 x2 11.76471
the determinant is even smaller
D 0.77(1) 1(0.705882) 0.064
9.4 (a) The determinant can be computed as:
A1 1 1 1(0) 1(1) 1
1 0
A2 2 1 2(0) 1(3) 3
3 0
A3 2 1 2(1) 1(3) 1
3 1
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
3
D 0(1) 2(3) 5(1) 1
(b) Cramer’s rule
1 2 5
1 1 1
2 1 0 2
x1 2
D 1
0 1 5
2 1 1
3 2 0 8
x2 8
D 1
0 2 1
2 1 1
3 1 2 3
x3 3
D 1
(c) Pivoting is necessary, so switch the first and third rows,
3x1 x2 2
2 x1 x2 x3 1
2 x2 5 x3 1
Multiply pivot row 1 by 2/3 and subtract the result from the second row to eliminate the a21 term. Note that
because a31 = 0, it does not have to be eliminated
3x1 x2 2
0.33333x2 x3 0.33333
2 x2 5 x3 1
Pivoting is necessary so switch the second and third row,
3x1 x2 2
2 x2 5 x3 1
0.33333x2 x3 0.33333
Multiply pivot row 2 by 0.33333/2 and subtract the result from the third row to eliminate the a32 term.
3x1 x2 2
2 x2 5 x3 1
0.16667 x3 0.5
The solution can then be obtained by back substitution
0.5
x3 3
0.16667
1 5(3)
x2 8
2
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
4
2 0(3) 1(8)
x1 2
3
The results can be checked by substituting them back into the original equations:
2(8) 5(3) 1
2(2) 8 3 1
3(2) 8 2
9.5 (a) The equations can be expressed in a format that is compatible with graphing x2 versus x1:
x2 0.5 x1 9.5
x2 0.51x1 9.4
The resulting plot indicates that the intersection of the lines is difficult to detect:
22
20
18
16
14
12
10
5 10 15 20
Only when the plot is zoomed is it at all possible to discern that solution seems to lie at about x1 = 14.5 and
x2 = 10.
14.7
14.65
14.6
14.55
14.5
14.45
14.4
14.35
14.3
9.75 10 10.25
(b) The determinant can be computed as
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
5
0.5 1
0.5(2) (1)(1.02) 0.02
1.02 2
which is close to zero.
(c) Because the lines have very similar slopes and the determinant is so small, you would expect that the
system would be ill-conditioned
(d) Multiply the first equation by 1.02/0.5 and subtract the result from the second equation to eliminate the
x1 term from the second equation,
0.5 x1 x2 9.5
0.04 x2 0.58
The second equation can be solved for
0.58
x2 14.5
0.04
This result can be substituted into the first equation which can be solved for
9.5 14.5
x1 10
0.5
(e) Multiply the first equation by 1.02/0.52 and subtract the result from the second equation to eliminate the
x1 term from the second equation,
0.52 x1 x2 9.5
0.03846 x2 0.16538
The second equation can be solved for
0.16538
x2 4.3
0.03846
This result can be substituted into the first equation which can be solved for
9.5 4.3
x1 10
0.52
Interpretation: The fact that a slight change in one of the coefficients results in a radically different solution
illustrates that this system is very ill-conditioned.
9.6 (a) Multiply the first equation by –3/10 and subtract the result from the second equation to eliminate the
x1 term from the second equation. Then, multiply the first equation by 1/10 and subtract the result from the
third equation to eliminate the x1 term from the third equation.
10 x1 2 x2 x3 27
4.4 x2 1.7 x3 53.4
0.8 x2 6.1x3 24.2
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
6
Multiply the second equation by 0.8/(–4.4) and subtract the result from the third equation to eliminate the x2
term from the third equation,
10 x1 2 x2 x3 27
4.4 x2 1.7 x3 53.4
6.409091x3 33.9091
Back substitution can then be used to determine the unknowns
33.9091
x3 5.29078
6.409091
(53.4 1.7(5.29078))
x2 10.0922
4.4
(27 5.29078 2(10.0922))
x1 0.152482
10
(b) Check:
10(0.152482) 2(10.0922) (5.29078) 27
3(0.152482) 5(10.0922) 2(5.29078) 61.5
0.152482 10.0922 5( 5.29078) 21.5
9.7 (a) Pivoting is necessary, so switch the first and third rows,
8 x1 x2 2 x3 20
3 x1 x2 7 x3 34
2 x1 6 x2 x3 38
Multiply the first equation by –3/(–8) and subtract the result from the second equation to eliminate the a21
term from the second equation. Then, multiply the first equation by 2/(–8) and subtract the result from the
third equation to eliminate the a31 term from the third equation.
8 x1 x2 2 x3 20
1.375 x2 7.75 x3 26.5
5.75 x2 1.5 x3 43
Pivoting is necessary so switch the second and third row,
8 x1 x2 2 x3 20
5.75 x2 1.5 x3 43
1.375 x2 7.75 x3 26.5
Multiply pivot row 2 by –1.375/(–5.75) and subtract the result from the third row to eliminate the a32 term.
8 x1 x2 2 x3 20
5.75 x2 1.5 x3 43
8.108696x3 16.21739
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
7
At this point, the determinant can be computed as
D 8 5.75 8.108696 ( 1) 2 373
The solution can then be obtained by back substitution
16.21739
x3 2
8.108696
43 1.5(2)
x2 8
5.75
20 2(2) 1(8)
x1 4
8
(b) Check:
2(4) 6(8) (2) 38
3(4) (8) 7(2) 34
8(4) (8) 2(2) 20
9.8 Multiply the first equation by –0.4/0.8 and subtract the result from the second equation to eliminate the
x1 term from the second equation.
0.8 0.4 x1 41
0.6 0.4 x2 45.5
0.4 0.8 x3 105
Multiply pivot row 2 by –0.4/0.6 and subtract the result from the third row to eliminate the x2 term.
0.8 0.4 x1 41
0.6 0.4 x2 45.5
0.533333 x3 135.3333
The solution can then be obtained by back substitution
135.3333
x3 253.75
0.533333
45.5 (0.4)253.75
x2 245
0.6
41 (0.4)245
x1 173.75
0.8
(b) Check:
0.8(173.75) 0.4(245) 41
0.4(173.75) 0.8(245) 0.4(253.75) 25
0.4(245) 0.8(253.75) 105
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
8
9.9 Mass balances can be written for each of the reactors as
200 Q13c1 Q12 c1 Q21c2 0
Q12 c1 Q21c2 Q23c2 0
500 Q13c1 Q23c2 Q33c3 0
Values for the flows can be substituted and the system of equations can be written in matrix form as
130 30 0 c1 200
90 90
0 c2 0
40 60 120 c3 500
The solution can then be developed using MATLAB,
>> A=[130 -30 0;-90 90 0;-40 -60 120];
>> B=[200;0;500];
>> C=A\B
C =
2.0000
2.0000
5.8333
9.10 Let xi = the volume taken from pit i. Therefore, the following system of equations must hold
0.52 x1 0.20 x2 0.25 x3 4800
0.30 x1 0.50 x2 0.20 x3 5800
0.18 x1 0.30 x2 0.55 x3 5700
MATLAB can be used to solve this system of equations for
>> A=[0.55 0.2 0.25;0.3 0.5 0.2;0.18 0.3 0.55];
>> b=[4800;5800;5700];
>> x=A\b
x =
1.0e+003 *
3.7263
7.2991
5.1628
Therefore, we take x1 = 3726.3, x2 = 7299.1, and x3 = 5162.8 m3 from pits 1, 2 and 3 respectively.
9.11 Let ci = component i. Therefore, the following system of equations must hold
15c1 17c2 19c3 2120
0.25c1 0.33c2 0.42c3 43.4
1.0c1 1.2c2 1.6c3 164
The solution can be developed with MATLAB:
A=[15 17 19;0.25 0.33 0.42;1 1.2 1.6];
b=[2120;43.4;164];
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
9
c=A\b
c =
20.0000
40.0000
60.0000
Therefore, c1 = 20, c2 = 40, and c3 = 60.
9.12 Centered differences (recall Chap. 4) can be substituted for the derivatives to give
ci 1 2ci ci 1 c c
0D U i 1 i 1 kci
x 2 2x
collecting terms yields
– (D + 0.5Ux)ci–1 + (2D + kx2)ci – (D – 0.5Ux)ci+1 = 0
Assuming x = 1 and substituting the parameters gives
2.5ci 1 4.2ci 1.5ci 1 0
For the first interior node (i = 1),
4.2c1 1.5c2 200
For the last interior node (i = 9)
2.5c8 4.2c9 15
These and the equations for the other interior nodes can be assembled in matrix form as
4.2 1.5 0 0 0 0 0 0 0 c1 200
2.5 4.2 1.5
0 0 0 0 0 0 c2 0
0 2.5 4.2 1.5 0 0 0 0 0 c3 0
0 0 2.5 4.2 1.5 0 0 0 0 c4 0
0
0 0 2.5 4.2 1.5 0 0 0 c5 0
0 0 0 0 2.5 4.2 1.5 0 0 c6 0
0 0 0 0 0 2.5 4.2 1.5 0 c7 0
0 0 0 0 0 0 2.5 4.2 1.5 c8 0
0 0 0 0 0 0 0 2.5 4.2 c9 15
The following script generates and plots the solution:
A=[4.2 -1.5 0 0 0 0 0 0 0
-2.5 4.2 -1.5 0 0 0 0 0 0
0 -2.5 4.2 -1.5 0 0 0 0 0
0 0 -2.5 4.2 -1.5 0 0 0 0
0 0 0 -2.5 4.2 -1.5 0 0 0
0 0 0 0 -2.5 4.2 -1.5 0 0
0 0 0 0 0 -2.5 4.2 -1.5 0
0 0 0 0 0 0 -2.5 4.2 -1.5
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
10
0 0 0 0 0 0 0 -2.5 4.2];
b=[200 0 0 0 0 0 0 0 15]';
c=A\b;
c=[80 c' 10];
x=0:1:10;
plot(x,c)
80
60
40
20
0
0 2 4 6 8 10
9.13 For the first stage, the mass balance can be written as
F1 yin F2 x2 F2 x1 F1 x1
Substituting x = Ky and rearranging gives
F F
1 2 K y1 2 Ky2 yin
F1 F1
Using a similar approach, the equation for the last stage is
F F
y4 1 2 K y5 2 xin
F1 F1
For interior stages,
F F
yi-1 1 2 K yi 2 Kyi 1 0
F1 F1
These equations can be used to develop the following system,
11 10 0 0 0 y1 0.1
1 11 10 0 0 y2 0
0 1 11 10 0 y3 0
0 0 1 11 10 y4 0
0 0 0 1 11 y5 0
The solution can be developed in a number of ways. For example, using MATLAB,
>> format short g
>> A=[11 -10 0 0 0;-1 11 -10 0 0;0 -1 11 -10 0;0 0 -1 11 -10;0 0 0 -1 11];
>> B=[0.1;0;0;0;0];
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
11
>> Y=A\B
Y =
0.0099999
0.0009999
9.99e-005
9.9e-006
9e-007
Note that the corresponding values of X can be computed as
>> X=5*Y
X =
0.05
0.0049995
0.0004995
4.95e-005
4.5e-006
Therefore, yout = 0.0000009 and xout = 0.05.
9.14 Assuming a unit flow for Q1, the simultaneous equations can be written in matrix form as
2 1 2 0 0 0 Q2 0
0 0 2 1 2 0 Q3 0
0 0 0 0 2 3 Q4 0
1 1 0 0 0 0 Q5 1
0 1 1 1 0 0 Q6 0
0 0 0 1 1 1 Q 0
7
These equations can then be solved with MATLAB,
>> A=[-2 1 2 0 0 0;
0 0 -2 1 2 0;
0 0 0 0 -2 3;
1 1 0 0 0 0;
0 1 -1 -1 0 0;
0 0 0 1 -1 -1];
>> B=[0 0 0 1 0 0 ]';
>> Q=A\B
Q =
0.5059
0.4941
0.2588
0.2353
0.1412
0.0941
9.15 The solution can be generated with MATLAB,
>> A=[1 0 0 0 0 0 0 0 1 0;
0 0 1 0 0 0 0 1 0 0;
0 1 0 3/5 0 0 0 0 0 0;
-1 0 0 -4/5 0 0 0 0 0 0;
0 -1 0 0 0 0 3/5 0 0 0;
0 0 0 0 -1 0 -4/5 0 0 0;
0 0 -1 -3/5 0 1 0 0 0 0;
0 0 0 4/5 1 0 0 0 0 0;
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
12
0 0 0 0 0 -1 -3/5 0 0 0;
0 0 0 0 0 0 4/5 0 0 1];
>> B=[0 0 -74 0 0 24 0 0 0 0]';
>> x=A\B
x =
37.3333
-46.0000
74.0000
-46.6667
37.3333
46.0000
-76.6667
-74.0000
-37.3333
61.3333
Therefore, in kN
AB = 37.3333 BC = 46 AD = 74 BD = 46.6667 CD = 37.3333
DE = 46 CE = 76.6667 Ax = 74 Ay = 37.33333 Ey = 61.3333
9.16
function x=pentasol(A,b)
% pentasol: pentadiagonal system solver banded system
% x=pentasol(A,b):
% Solve a pentadiagonal system Ax=b
% input:
% A = pentadiagonal matrix
% b = right hand side vector
% output:
% x = solution vector
% Error checks
[m,n]=size(A);
if m~=n,error('Matrix must be square');end
if length(b)~=m,error('Matrix and vector must have the same number of
rows');end
x=zeros(n,1);
% Extract bands
d=[0;0;diag(A,-2)];
e=[0;diag(A,-1)];
f=diag(A);
g=diag(A,1);
h=diag(A,2);
delta=zeros(n,1);
epsilon=zeros(n-1,1);
gamma=zeros(n-2,1);
alpha=zeros(n,1);
c=zeros(n,1);
z=zeros(n,1);
% Decomposition
delta(1)=f(1);
epsilon(1)=g(1)/delta(1);
gamma(1)=h(1)/delta(1);
alpha(2)=e(2);
delta(2)=f(2)-alpha(2)*epsilon(1);
epsilon(2)=(g(2)-alpha(2)*gamma(1))/delta(2);
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
13
gamma(2)=h(2)/delta(2);
for k=3:n-2
alpha(k)=e(k)-d(k)*epsilon(k-2);
delta(k)=f(k)-d(k)*gamma(k-2)-alpha(k)*epsilon(k-1);
epsilon(k)=(g(k)-alpha(k)*gamma(k-1))/delta(k);
gamma(k)=h(k)/delta(k);
end
alpha(n-1)=e(n-1)-d(n-1)*epsilon(n-3);
delta(n-1)=f(n-1)-d(n-1)*gamma(n-3)-alpha(n-1)*epsilon(n-2);
epsilon(n-1)=(g(n-1)-alpha(n-1)*gamma(n-2))/delta(n-1);
alpha(n)=e(n)-d(n)*epsilon(n-2);
delta(n)=f(n)-d(n)*gamma(n-2)-alpha(n)*epsilon(n-1);
% Forward substitution
c(1)=b(1)/delta(1);
c(2)=(b(2)-alpha(2)*c(1))/delta(2);
for k=3:n
c(k)=(b(k)-d(k)*c(k-2)-alpha(k)*c(k-1))/delta(k);
end
% Back substitution
x(n)=c(n);
x(n-1)=c(n-1)-epsilon(n-1)*x(n);
for k=n-2:-1:1
x(k)=c(k)-epsilon(k)*x(k+1)-gamma(k)*x(k+2);
end
Test of function:
>> A=[8 -2 -1 0 0
-2 9 -4 -1 0
-1 3 7 -1 -2
0 -4 -2 12 -5
0 0 -7 -3 15];
>> b=[5 2 1 1 5]';
>> x=pentasol(A,b)
x =
0.7993
0.5721
0.2503
0.5491
0.5599
9.17 Here is the M-file function based on Fig. 9.5 to implement Gauss elimination with partial pivoting
function [x, D] = GaussPivotNew(A, b, tol)
% GaussPivotNew: Gauss elimination pivoting
% [x, D] = GaussPivotNew(A,b,tol): Gauss elimination with pivoting.
% input:
% A = coefficient matrix
% b = right hand side vector
% tol = tolerance for detecting "near zero"
% output:
% x = solution vector
% D = determinant
[m,n]=size(A);
if m~=n, error('Matrix A must be square'); end
nb=n+1;
Aug=[A b];
npiv=0;
% forward elimination
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
14
for k = 1:n-1
% partial pivoting
[big,i]=max(abs(Aug(k:n,k)));
ipr=i+k-1;
if ipr~=k
npiv=npiv+1;
Aug([k,ipr],:)=Aug([ipr,k],:);
end
absakk=abs(Aug(k,k));
if abs(Aug(k,k))<=tol
D=0;
error('Singular or near singular system')
end
for i = k+1:n
factor=Aug(i,k)/Aug(k,k);
Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);
end
end
for i = 1:n
if abs(Aug(i,i))<=tol
D=0;
error('Singular or near singular system')
end
end
% back substitution
x=zeros(n,1);
x(n)=Aug(n,nb)/Aug(n,n);
for i = n-1:-1:1
x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
end
D=(-1)^npiv;
for i=1:n
D=D*Aug(i,i);
end
Here is a script to solve Prob. 9.5 for the two cases of tol:
clear; clc; format short g
A=[0.5 -1;1.02 -2];
b=[-9.5;-18.8];
disp('Solution and determinant calculated with GaussPivotNew:')
[x, D] = GaussPivotNew(A,b,1e-5)
disp('Determinant calculated with det:')
D=det(A)
The resulting output is
Solution and determinant calculated with GaussPivotNew:
x =
10
14.5
D =
0.02
Determinant calculated with det:
D =
0.02
PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual
may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the
publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their
individual course preparation. If you are a student using this Manual, you are using it without permission.
