Chapra Numerical Analysis 11th Chapter Solution
Chapra Numerical Analysis 11th Chapter Solution
CHAPTER 11
11.1 First, the decomposition is implemented as
e2 = 0.4/0.8 = 0.5
f2 = 0.8 0.5)(0.4) = 0.6
e3 = 0.4/0.6 = 0.66667
f3 = 0.8 0.66667)(0.4) = 0.53333
Transformed system is
0.4
0.8
0.6
0.5
0
0.66667
0
0.4
0.53333
which is decomposed as
0
1
[L] 0.5
1
0
0.66667
0
0
1
0
0.8 0.4
[U ] 0
0.6
0.4
0
0
0.53333
0.49
1
[L]
0.645
1
0.717 1
2.04 1
1.550 1
[U ]
1.395 1
1.323
1
[ L]{D} 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.
2
Solving this gives
0.490196
{D}
0.316296
0.226775
Back substitution, [U]{X} = {D}, can then be implemented to give to first column of the
inverse
0.755841
{X } 0.541916
0.349667
0.171406
0
[ L]{D} 1
0
0
which leads to
0.541916
{X } 1.105509
0.713322
0.349667
0
[ L]{D} 0
1
0
which leads to
0.349667
{X } 0.713322
1.105509
0.541916
0
[ L]{D} 0
0
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.
which leads to
0.171406
{X } 0.349667
0.541916
0.755841
0.755841
[ A] 0.541916
0.349667
0.171406
1
0.541916
1.105509
0.713322
0.349667
0.349667
0.713322
1.105509
0.541916
0.171406
0.349667
0.541916
0.755841
2.01475 0.02875
0.01036 2.014534
which is decomposed as
0.01036
1
[L]
0.01036
1
0.01036 1
2.01475 0.02875
2.014534 0.02875
[U ]
2.014534 0.02875
2.014534
4
r4 = 2.087505
Back substitution
x4 = 1.036222
x3 = 0.01096
x2 = 0.021586
x1 = 2.072441
11.4 We can use MATLAB to verify the results of Example 11.2,
>> L=[2.4495 0 0;6.1237 4.1833 0;22.454 20.916 6.1106]
L =
2.4495
6.1237
22.4540
0
4.1833
20.9160
0
0
6.1106
15.0000
54.9997
224.9995
55.0011
224.9995
979.0006
>> L*L'
ans =
6.0001
15.0000
55.0011
11.5
l11 8 2.828427
l 21
20
7.071068
2.828427
l 22 80 7.071068 2 5.477226
l31
15
5.303301
2.828427
l32
50 7.071068(5.303301)
2.282177
5.477226
l 33 60 5.303301 2 2.282177
5.163978
2.828427
11.6
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
l11 6 2.44949
l 21
15
6.123724
2.44949
l 22 55 6.123724 2 4.1833
l31
55
22.45366
2.44949
l32
225 6.123724(22.45366)
20.9165
4.1833
2.44949
The solution can then be generated by first using forward substitution to modify the righthand-side vector,
[ L]{D} {B}
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.
x1
41 0.4 x2 41 0.4(0)
51.25
0.8
0.8
x2
56.875
0.8
0.8
x3
159.6875
0.8
0.8
Second iteration:
x1
41 0.4(56.875)
79.6875
0.8
x2
25 0.4(79.6875) 0.4(159.6875)
150.9375
0.8
x3
105 0.4(150.9375)
206.7188
0.8
a,1
79 .6875 51 .25
100 % 35 .69 %
79 .6875
a,2
a ,3
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
1
unknown
x1
x2
x3
x1
x2
x3
x1
x2
x3
x1
x2
value
51.25
56.875
159.6875
79.6875
150.9375
206.7188
126.7188
197.9688
230.2344
150.2344
221.4844
a
100.00%
100.00%
100.00%
35.69%
62.32%
22.75%
37.11%
23.76%
10.21%
15.65%
10.62%
maximum a
100.00%
62.32%
37.11%
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.
x3
x1
x2
x3
x1
x2
x3
241.9922
161.9922
233.2422
247.8711
167.8711
239.1211
250.8105
4.86%
7.26%
5.04%
2.37%
3.50%
2.46%
1.17%
15.65%
7.26%
3.50%
Thus, after 6 iterations, the maximum error is 3.5% and we arrive at the result: x1 = 167.8711,
x2 = 239.1211 and x3 = 250.8105.
(b) The same computation can be developed with relaxation where = 1.2.
First iteration:
x1
41 0.4 x 2 41 0.4(0)
51.25
0.8
0.8
x2
62
0.8
0.8
x3
168.45
0.8
0.8
x1
41 0.4(74.4)
88.45
0.8
x2
25 0.4(93.84) 0.4(202.14)
179.24
0.8
x3
105 0.4(200.208)
231.354
0.8
a ,1
93 .84 61 .5
100 % 34 .46 %
93 .84
a,2
200 .208 74 .4
100 % 62 .84 %
200 .208
a ,3
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
1
unknown
x1
x2
x3
x1
x2
x3
x1
x2
x3
x1
x2
x3
value
51.25
62
168.45
88.45
179.24
231.354
151.354
231.2768
249.99528
169.99528
243.23898
253.44433
relaxation
61.5
74.4
202.14
93.84
200.208
237.1968
162.8568
237.49056
252.55498
171.42298
244.38866
253.6222
a
100.00%
100.00%
100.00%
34.46%
62.84%
14.78%
42.38%
15.70%
6.08%
5.00%
2.82%
0.42%
maximum a
100.000%
62.839%
42.379%
4.997%
Thus, relaxation speeds up convergence. After 6 iterations, the maximum error is 4.997% and
we arrive at the result: x1 = 171.423, x2 = 244.389 and x3 = 253.622.
11.8 The first iteration can be implemented as
c1
253.3333
15
15
c2
108.8889
18
18
c3
289.3519
12
12
Second iteration:
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.
c1
294.4012
15
15
c2
212.1842
18
18
c3
311.6491
12
12
a ,1
a,2
a ,3
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
1
unknown
c1
c2
c3
c1
c2
c3
c1
c2
c3
c1
c2
c3
value
253.3333
108.8889
289.3519
294.4012
212.1842
311.6491
316.5468
223.3075
319.9579
319.3254
226.5402
321.1535
a
100.00%
100.00%
100.00%
13.95%
48.68%
7.15%
7.00%
4.98%
2.60%
0.87%
1.43%
0.37%
maximum a
100.00%
48.68%
7.00%
1.43%
Thus, after 4 iterations, the maximum error is 1.43% and we arrive at the result: c1 =
319.3254, c2 = 226.5402 and c3 = 321.1535.
11.9 The first iteration can be implemented as
c1
253.3333
15
15
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
c2
66.6667
18
18
c3
195.8333
12
12
Second iteration:
c1
279.7222
15
15
c2
174.1667
18
18
c3
285.8333
12
12
a ,1
a,2
a ,3
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
1
unknown
c1
c2
c3
c1
c2
c3
c1
c2
c3
c1
c2
c3
value
253.3333
66.66667
195.8333
279.7222
174.1667
285.8333
307.2222
208.5648
303.588
315.2855
219.0664
315.6211
a
100.00%
100.00%
100.00%
9.43%
61.72%
31.49%
8.95%
16.49%
5.85%
2.56%
4.79%
3.81%
maximum a
100.00%
61.72%
16.49%
4.79%
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
Thus, after 4 iterations, the maximum error is 4.79% and we arrive at the result: c1 =
315.5402, c2 = 219.0664 and c3 = 315.6211.
11.10 The first iteration can be implemented as
x1
27 2 x 2 x3 27 2(0) 0
2.7
10
10
x2
8.9
6
6
x3
6.62
5
5
Second iteration:
x1
27 2(8.9) 6.62
0.258
10
x2
x3
a,1
0.258 2.7
100 % 947 %
0.258
a,2
7.914333 8.9
100 % 12 .45 %
7.914333
a ,3
5.934467 (6.62 )
100 % 11 .55 %
5.934467
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
1
unknown
x1
x2
x3
x1
x2
x3
value
2.7
8.9
-6.62
0.258
7.914333
-5.93447
a
100.00%
100.00%
100.00%
946.51%
12.45%
11.55%
maximum a
100%
946%
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
x1
x2
x3
x1
x2
x3
x1
x2
x3
0.523687
8.010001
-6.00674
0.497326
7.999091
-5.99928
0.500253
8.000112
-6.00007
50.73%
1.19%
1.20%
5.30%
0.14%
0.12%
0.59%
0.01%
0.01%
50.73%
5.30%
0.59%
Thus, after 5 iterations, the maximum error is 0.59% and we arrive at the result: x1 =
0.500253, x2 = 8.000112 and x3 = 6.00007.
11.11 The equations should first be rearranged so that they are diagonally dominant,
6 x1 x 2 x3 3
6 x1 9 x 2 x3 40
3x1 x 2 12 x3 50
Each can be solved for the unknown on the diagonal as
x1
3 x 2 x3
6
x2
40 6 x1 x3
9
x3
50 3x1 x 2
12
300
0.5
6
x2
40 6(0.5) 0
4.11111
9
x3
50 3(0.5) 4.11111
3.949074
12
Second iteration:
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
x1
3 4.11111 3.949074
1.843364
6
x2
40 6(1.843364 ) 3.949074
2.776749
9
x3
50 3(1.843364 ) 2.776749
4.396112
12
a ,1
1.843364 0.5
100 % 72 .88 %
1.843364
a,2
2.776749 4.11111
100 % 48 .05 %
2.776749
a ,3
4.396112 3.949074
100 % 10 .17 %
4.396112
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
1
unknown
x1
x2
x3
x1
x2
x3
x1
x2
x3
x1
x2
x3
value
0.5
4.111111
3.949074
1.843364
2.776749
4.396112
1.695477
2.82567
4.355063
1.696789
2.829356
4.355084
a
100.00%
100.00%
100.00%
72.88%
48.05%
10.17%
8.72%
1.73%
0.94%
0.08%
0.13%
0.00%
maximum a
100.00%
72.88%
8.72%
0.13%
Thus, after 4 iterations, the maximum error is 0.13% and we arrive at the result: x1 =
1.696789, x2 = 2.829356 and x3 = 4.355084.
(b) First iteration: To start, assume x1 = x2 = x3 = 0
x1new
300
0.5
6
Apply relaxation
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
x1 0.95(0.5) (1 0.95 )0 0.475
40 6(0.475) 0
4.12778
9
x2new
50 3(0.475) 3.92139
3.95863
12
x3new
Note that error estimates are not made on the first iteration, because all errors will be 100%.
Second iteration:
x1new
3 3.92139 3.76070
1.78035
6
a ,1
1.71508 0.475
100 % 72 .3%
1.71508
Because this error exceeds the stopping criterion, it will not be necessary to compute error
estimates for the remainder of this iteration.
x2new
40 6(1.71508) 3.76070
2.88320
9
x3new
50 3(1.71508) 2.93511
4.35084
12
The computations can be continued for one more iteration. The entire calculation is
summarized in the following table.
iteration
1
2
3
x1
0.50000
1.78035
1.70941
x1r
0.47500
1.71508
1.70969
a1
100.0%
72.3%
0.3%
x2
4.12778
2.88320
2.82450
x2r
3.92139
2.93511
2.83003
a2
100.0%
33.6%
3.7%
x3
3.95863
4.35084
4.35825
x3r
3.76070
4.32134
4.35641
a3
100.0%
13.0%
0.8%
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.
15
After 3 iterations, the approximate errors fall below the stopping criterion with the final
result: x1 = 1.70969, x2 = 2.82450 and x3 = 4.35641. Note that the exact solution is x1 =
1.69737, x2 = 2.82895 and x3 = 4.35526
11.12 The equations must first be rearranged so that they are diagonally dominant
8 x1 x 2 2 x3 20
2 x1 6 x 2 x3 38
3x1 x 2 7 x3 34
x1
20 x2 2 x3 20 0 2(0)
2.5
8
8
x2
38 2 x1 x3 38 2(2.5) 0
7.166667
6
6
x3
2.761905
7
7
Second iteration:
x1
20 7.166667 2(2.761905)
4.08631
8
x2
38 2 x1 x3 38 2(4.08631) (2.761905)
8.155754
6
6
x3
1.94076
7
7
a ,1
4.08631 2.5
100 % 38 .82 %
4.08631
a,2
8.155754 7.166667
100 % 12 .13 %
8.155754
a ,3
1.94076 (2.761905 )
100 % 42 .31 %
1.94076
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.
16
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
0
unknown
x1
x2
x3
x1
x2
x3
x1
x2
x3
x1
x2
x3
value
0
0
0
2.5
7.166667
-2.7619
4.08631
8.155754
-1.94076
4.004659
7.99168
-1.99919
100.00%
100.00%
100.00%
38.82%
12.13%
42.31%
2.04%
2.05%
2.92%
maximum a
100.00%
42.31%
2.92%
Thus, after 3 iterations, the maximum error is 2.92% and we arrive at the result: x1 =
4.004659, x2 = 7.99168 and x3 = 1.99919.
(b) The same computation can be developed with relaxation where = 1.2.
First iteration:
x1
20 x2 2 x3 20 0 2(0)
2.5
8
8
x2
38 2 x1 x3 38 2(3) 0
7.333333
6
6
x3
2.3142857
7
7
x1
20 x 2 2 x3 20 8.8 2(2.7771429)
4.2942857
8
8
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.
17
x2
38 2 x1 x3 38 2(4.5531429) 2.7771429
8.3139048
6
6
x3
1.7319837
7
7
a ,1
4.5531429 3
100 % 34 .11 %
4.5531429
a,2
8.2166857 8.8
100 % 7.1%
8.2166857
a ,3
1.5229518 (2.7771429 )
100 % 82 .35 %
1.5229518
The remainder of the calculation proceeds until all the errors fall below the stopping criterion
of 5%. The entire computation can be summarized as
iteration
1
unknown
x1
x2
x3
x1
x2
x3
x1
x2
x3
x1
x2
x3
x1
x2
x3
x1
x2
x3
value
2.5
7.3333333
-2.314286
4.2942857
8.3139048
-1.731984
3.9078237
7.8467453
-2.12728
4.0336312
8.0695595
-1.945323
3.9873047
7.9700747
-2.022594
4.0048286
8.0124354
-1.990866
relaxation
3
8.8
-2.777143
4.5531429
8.2166857
-1.522952
3.7787598
7.7727572
-2.248146
4.0846055
8.12892
-1.884759
3.9678445
7.9383056
-2.050162
4.0122254
8.0272613
-1.979007
a
100.00%
100.00%
100.00%
34.11%
7.10%
82.35%
20.49%
5.71%
32.26%
7.49%
4.38%
19.28%
2.94%
2.40%
8.07%
1.11%
1.11%
3.60%
maximum a
100.000%
82.353%
32.257%
19.280%
8.068%
3.595%
Thus, relaxation actually seems to retard convergence. After 6 iterations, the maximum error
is 3.595% and we arrive at the result: x1 = 4.0122254, x2 = 8.0272613 and x3 = 1.979007.
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.
18
11.13 As shown below, for slopes of 1 and 1 the Gauss-Seidel technique will neither converge
nor diverge but will oscillate interminably.
x2
v
x1
u
11.14 As ordered, none of the sets will converge. However, if Set 1 and 2 are reordered so that
they are diagonally dominant, they will converge on the solution of (1, 1, 1).
Set 1:
9x + 3y + z = 13
2x + 5y z = 6
6x
+ 8z = 2
Set 2:
4x + 2y z = 4
x + 5y z = 5
x + y + 6z = 8
At face value, because it is not strictly diagonally dominant, Set 2 would seem to be
divergent. However, since it is very close to being diagonally dominant, a solution can be
obtained.
The third set is not diagonally dominant and will diverge for most orderings. However, the
following arrangement will converge albeit at a very slow rate:
Set 3:
3x + 4y + 5z = 6
2y z = 1
2x + 2y 3z = 3
19
>> inv(A)
ans =
3.8750
-5.5000
2.1250
-5.5000
7.0000
-2.5000
2.1250
-2.5000
0.8750
>> cond(A,inf)
ans =
750.0000
(b) However, for the 44 system, the ill-conditioned nature of the matrix yields poor results:
>> A=[1 4 9 16;4 9 16 25;9 16 25 36;16 25 36 49];
>> B=[30 54 86 126]';
>> x=A\B
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 3.037487e-019.
x =
0.5496
2.3513
-0.3513
1.4504
>> cond(A,inf)
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 3.037487e-019.
> In cond at 48
ans =
3.2922e+018
Note that using other software such as Excel yields similar results. For example, the condition
number computed with Excel is 51017.
11.16 (a) As shown, there are 4 roots, one in each quadrant.
8
(0.618,3.236)
4
(1, 2)
0
-4
-2
0
-4
2
(1.618, 1.236)
(2, 4)
-8
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.
20
(b) It might be expected that if an initial guess was within a quadrant, the result would be the
root in the quadrant. However a sample of initial guesses spanning the range yield the
following roots:
6
3
0
-3
-6
(-2, -4)
(-0.618,3.236) (-0.618,3.236)
(1,2)
(-0.618,3.236) (-0.618,3.236) (-0.618,3.236)
(1,2)
(1,2)
(1.618, -1.236) (1.618, -1.236) (1.618, -1.236)
(-2, -4)
(-2, -4)
(1.618, -1.236) (1.618, -1.236)
(-2, -4)
(-2, -4)
(-2, -4)
(1.618, -1.236)
-6
-3
0
3
(-0.618,3.236)
(-0.618,3.236)
(1.618, -1.236)
(1.618, -1.236)
(-2, -4)
6
We have highlighted the guesses that converge to the roots in their quadrants. Although some
follow the pattern, others jump to roots that are far away. For example, the guess of (6, 0)
jumps to the root in the first quadrant.
This underscores the notion that root location techniques are highly sensitive to initial guesses
and that open methods like the Solver can locate roots that are not in the vicinity of the initial
guesses.
11.17 Define the quantity of transistors, resistors, and computer chips as x1, x2 and x3. The system
equations can then be defined as
4 x1 3x 2 2 x3 960
x1 3x 2 x3 510
2 x1 x 2 3x3 610
The solution can be implemented in Excel as shown below:
The following view shows the formulas that are employed to determine the inverse in cells
A7:C9 and the solution in cells D7:D9.
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.
21
Here is the same solution generated in MATLAB:
>> A=[4 3 2;1 3 1;2 1 3];
>> B=[960 510 610]';
>> x=A\B
x =
120
100
90
The digits of precision that could be lost due to ill-conditioning can be calculated as
>> c = log10(N)
c =
13.2048
Thus, about 13 digits could be suspect. A right-hand side vector can be developed
corresponding to a solution of ones:
>> b=[sum(A(1,:)); sum(A(2,:)); sum(A(3,:)); sum(A(4,:)); sum(A(5,:));
sum(A(6,:)); sum(A(7,:)); sum(A(8,:)); sum(A(9,:)); sum(A(10,:))]
b =
2.9290
2.0199
1.6032
1.3468
1.1682
1.0349
0.9307
0.8467
0.7773
0.7188
22
1.0000
1.0000
0.9999
1.0003
0.9995
1.0005
0.9997
1.0001
Thus, some of the results are accurate to only about 3 to 4 significant digits. Because
MATLAB represents numbers to 15 significant digits, this means that about 11 to 12 digits
are suspect.
11.19 First, the Vandermonde matrix can be set up
>> x1 = 4;x2=2;x3=7;x4=10;x5=3;x6=5;
>> A = [x1^5 x1^4 x1^3 x1^2 x1 1;x2^5 x2^4 x2^3 x2^2 x2 1;x3^5 x3^4
x3^3 x3^2 x3 1;x4^5 x4^4 x4^3 x4^2 x4 1;x5^5 x5^4 x5^3 x5^2 x5 1;x6^5
x6^4 x6^3 x6^2 x6 1]
A =
1024
32
16807
100000
243
3125
256
16
2401
10000
81
625
64
8
343
1000
27
125
16
4
49
100
9
25
4
2
7
10
3
5
1
1
1
1
1
1
The digits of precision that could be lost due to ill-conditioning can be calculated as
>> c = log10(N)
c =
7.1611
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.
23
Thus, about 7 digits might be suspect. A right-hand side vector can be developed
corresponding to a solution of ones:
>> b=[sum(A(1,:));sum(A(2,:));sum(A(3,:));sum(A(4,:));sum(A(5,:));
sum(A(6,:))]
b =
1365
63
19608
111111
364
3906
Some of the results are accurate to about 12 significant digits. Because MATLAB represents
numbers to about 15 significant digits, this means that about 3 digits are suspect. Thus, for
this case, the condition number tends to exaggerate the impact of ill-conditioning.
11.20 The flop counts for the tridiagonal algorithm in Fig. 11.2 can be determined as
Sub Decomp(e, f, g, n)
Dim k As Integer
For k = 2 To n
e(k) = e(k) / f(k - 1)
f(k) = f(k) - e(k) * g(k - 1)
Next k
End Sub
mult/div
add/subt
'(n 1)
'(n 1)
(n 1)
Sub Substitute(e, f, g, r, n, x)
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.
24
Dim k As Integer
For k = 2 To n
r(k) = r(k) - e(k) * r(k - 1)
Next k
x(n) = r(n) / f(n)
For k = n - 1 To 1 Step -1
x(k) = (r(k) - g(k) * x(k + 1)) / f(k)
Next k
End Sub
Sum =
'(n 1)
'
(n 1)
'2(n 1)
(n 1)
5(n-1) + 1
(3n 3)
The multiply/divides and add/subtracts can be summed to yield 8n 7 as opposed to n3/3 for
naive Gauss elimination. Therefore, a tridiagonal solver is well worth using.
1000000
100000
Tridiagonal
Naive Gauss
10000
1000
100
10
1
1
10
100
11.21 Here is a VBA macro to obtain a solution for a tridiagonal system using the Thomas
algorithm. It is set up to duplicate the results of Example 11.1.
Option Explicit
Sub TriDiag()
Dim i As Integer, n As Integer
Dim e(10) As Double, f(10) As Double, g(10) As Double
Dim r(10) As Double, x(10) As Double
n = 4
e(2) = -1: e(3) = -1: e(4) = -1
f(1) = 2.04: f(2) = 2.04: f(3) = 2.04: f(4) = 2.04
g(1) = -1: g(2) = -1: g(3) = -1
r(1) = 40.8: r(2) = 0.8: r(3) = 0.8: r(4) = 200.8
Call Thomas(e, f, g, r, n, x)
For i = 1 To n
MsgBox x(i)
Next i
End Sub
Sub Thomas(e, f, g, r, n, x)
Call Decomp(e, f, g, n)
Call Substitute(e, f, g, r, n, x)
End Sub
Sub Decomp(e, f, g, n)
Dim k As Integer
For k = 2 To n
e(k) = e(k) / f(k - 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.
25
f(k) = f(k) - e(k) * g(k - 1)
Next k
End Sub
Sub Substitute(e, f, g, r, n, x)
Dim k As Integer
For k = 2 To n
r(k) = r(k) - e(k) * r(k - 1)
Next k
x(n) = r(n) / f(n)
For k = n - 1 To 1 Step -1
x(k) = (r(k) - g(k) * x(k + 1)) / f(k)
Next k
End Sub
11.22 Here is a VBA macro to obtain a solution of a symmetric system with Cholesky
decomposition. It is set up to duplicate the results of Example 11.2.
Option Explicit
Sub TestChol()
Dim i As Integer, j As Integer
Dim n As Integer
Dim a(10, 10) As Double
n = 3
a(1, 1) = 6: a(1, 2) = 15: a(1, 3) = 55
a(2, 1) = 15: a(2, 2) = 55: a(2, 3) = 225
a(3, 1) = 55: a(3, 2) = 225: a(3, 3) = 979
Call Cholesky(a, n)
'output results to worksheet
Sheets("Sheet1").Select
Range("a3").Select
For i = 1 To n
For j = 1 To n
ActiveCell.Value = a(i, j)
ActiveCell.Offset(0, 1).Select
Next j
ActiveCell.Offset(1, -n).Select
Next i
Range("a3").Select
End Sub
Sub Cholesky(a, n)
Dim i As Integer, j As Integer, k As Integer
Dim sum As Double
For k = 1 To n
For i = 1 To k - 1
sum = 0
For j = 1 To i - 1
sum = sum + a(i, j) * a(k, j)
Next j
a(k, i) = (a(k, i) - sum) / a(i, i)
Next i
sum = 0
For j = 1 To k - 1
sum = sum + a(k, j) ^ 2
Next j
a(k, k) = Sqr(a(k, k) - sum)
Next k
End Sub
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.
26
11.23 Here is a VBA macro to obtain a solution of a linear diagonally-dominant system with the
Gauss-Seidel method. It is set up to duplicate the results of Example 11.3.
Option Explicit
Sub Gausseid()
Dim n As Integer, imax As Integer, i As Integer
Dim a(3, 3) As Double, b(3) As Double, x(3) As Double
Dim es As Double, lambda As Double
n = 3
a(1, 1) = 3: a(1, 2) = -0.1: a(1, 3) = -0.2
a(2, 1) = 0.1: a(2, 2) = 7: a(2, 3) = -0.3
a(3, 1) = 0.3: a(3, 2) = -0.2: a(3, 3) = 10
b(1) = 7.85: b(2) = -19.3: b(3) = 71.4
es = 0.1
imax = 20
lambda = 1#
Call Gseid(a, b, n, x, imax, es, lambda)
For i = 1 To n
MsgBox x(i)
Next i
End Sub
Sub Gseid(a, b, n, x, imax, es, lambda)
Dim i As Integer, j As Integer, iter As Integer, sentinel As Integer
Dim dummy As Double, sum As Double, ea As Double, old As Double
For i = 1 To n
dummy = a(i, i)
For j = 1 To n
a(i, j) = a(i, j) / dummy
Next j
b(i) = b(i) / dummy
Next i
For i = 1 To n
sum = b(i)
For j = 1 To n
If i <> j Then sum = sum - a(i, j) * x(j)
Next j
x(i) = sum
Next i
iter = 1
Do
sentinel = 1
For i = 1 To n
old = x(i)
sum = b(i)
For j = 1 To n
If i <> j Then sum = sum - a(i, j) * x(j)
Next j
x(i) = lambda * sum + (1# - lambda) * old
If sentinel = 1 And x(i) <> 0 Then
ea = Abs((x(i) - old) / x(i)) * 100
If ea > es Then sentinel = 0
End If
Next i
iter = iter + 1
If sentinel = 1 Or iter >= imax Then Exit Do
Loop
End Sub
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.