White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI
Подождите немного. Документ загружается.
68 CHAPTER 2. STEADY STATE DISCRETE MODELS
Since
E is an inverse of D and use the associative property,
E(D(E F)) = E0 = 0
(
ED)(E F) = 0
L (E F) = 0=
Proof that A
1
d is a solution of Ax = d:
Let
{ = D
1
g and again use the associative property
D(D
1
g) = (DD
1
)g = Lg = g=
Proofs of properties 3 and 4 are also a consequence of the associative prop-
erty.
2.2.7 Exercises
1. Consider the following algebraic system
1
{
1
+ 2{
2
+ 3{
3
= 1
1{
1
+ 1{
2
1{
3
= 2
2
{
1
+ 4{
2
+ 3{
3
= 3=
(a). Find the augmented matrix.
(b). By hand calculations with row operations and elementary matrices
find
H so that HD = X is upper triangular.
(c). Use this to find the solution, and verify your calculations using
MATLAB.
2. Use the M
ATLA B code heatgelm.m and experiment with the mesh sizes,
by using
q = 11> 21 and 41, in the heat conduction problem and verify that the
computed solution converges as the mesh goes to zero, that is, x
l
x(lk) goes
to zero as
k goes to zero
3. Prove property 3 of Theorem 2.2.1.
4. Prove property 4 of Theorem 2.2.1.
5. Prove that the solution of
D{ = g is unique if A
1
exists.
2.3 Cooling Fin and Tridiagonal Matrices
2.3.1 Introduction
In the thin wire problem we derived a tridiagonal matrix, which was gener-
ated from the finite di
erence approximation of the dierential equation. It
is very common to obtain either similar tridiagonal matrices or more compli-
cated matrices that have blocks of tridiagonal matrices. We will illustrate this
by a sequence of mo dels for a cooling fin. This section is concerned with a
very e
!cient version of the Gaussian elimination algorithm for the solution of
© 2004 by Chapman & Hall/CRC
2.3. COOLING FIN AND TRIDIAGONAL MATRICES 69
Figure 2.3.1: Thin Cooling Fin
tridiagonal algebraic systems. The full version of a Gaussian elimination al-
gorithm for
q unknowns requires order q
3
@3 operations and order q
2
storage
locations. By taking advantage of the number of zeros and their location, the
Gaussian elimination algorithm for tridiagonal systems can be reduced to order
5
q operations and order 8q storage locations!
2.3.2 Applied Area
Consider a hot mass, which must be cooled by transferring heat from the mass
to a cooler surrounding region. Examples include computer chips, electrical
amplifiers, a transformer on a power line, or a gasoline engine. One way to do
this is to attach cooling fins to this mass so that the surface area that transmits
the heat will be larger. We wish to be able to model heat flow so that one can
determine whether or not a particular configuration will su
!ciently cool the
mass.
In order to start the modeling process, we will make some assumptions that
will simplify the model. Later we will return to this model and reconsider some
of these assumptions. First, assume no time dependence and the temperature
is approximated by a function of only the distance from the mass to be cooled.
Thus, there is d i
usion in only one direction. This is depicted in Figure 2.3.1
where
{ is the direction perpendicular to the hot mass.
Second, assume the heat lost through the surface of the fin is similar to
Newton’s law of cooling so that for a slice of the lateral surface
heat loss through a slice = (area)(time interval)
f(x
vxu
x)
=
k(2Z + 2W ) w f(x
vxu
x)=
Here x
vxu
is the surrounding temp erature, and the f reflects the ability of the
fin’s surface to transmit heat to the surrounding region. If
f is n ear zero, then
© 2004 by Chapman & Hall/CRC
70 CHAPTER 2. STEADY STATE DISCRETE MODELS
little heat is lost. If
f is large, then a larger amount of heat is lost through the
lateral surface.
Third, assume heat di
uses in the { direction according to Fourier’s heat
law where N is the thermal conductivity. For interior volume elements with
{ ? O = 1,
0
(heat through lateral surface )
+(heat di
using through front)
(heat diusing through back)
=
k (2Z + 2W ) w f(x
vxu
x({))
+
W Z w Nx
{
({ + k@2)
W Z w Nx
{
({ k@2)= (2.3.1)
For the tip of the fin with
{ = O, we use Nx
{
(O) = f(x
vxu
x(O)) and
0
(heat through lateral surface of tip)
+(heat di
using through front)
(heat diusing through back)
= (
k@2)(2Z + 2W ) w f(x
vxu
x(O))
+
W Z w f(x
vxu
x(O))
W Z w Nx
{
(O k@2)= (2.3.2)
Note, the volume element near the tip of the fin is one half of the volume of the
interior elements.
These are only approximations because the temperature changes continu-
ously with space. In order to make these approximations in (2.3.1) and (2.3.2)
more accurate, we divide by
k w W Z and let k go to zero
0 = (2
Z + 2W )@(W Z ) f(x
vxu
x) + (Nx
{
)
{
= (2.3.3)
Let
F ((2Z + 2W )@(W Z )) f and i Fx
vxu
. The continuous model is given
by the following di
erential equation and two boundary conditions.
(Nx
{
)
{
+ Fx = i> (2.3.4)
x(0) = given and (2.3.5)
Nx
{
(O) = f(x
vxu
x(O))= (2.3.6)
The boundary condition in (2.3.6) is often called a derivative or flux or Robin
boundary condition.. If
f = 0, then no heat is allowed to pass through the
right boundary, and this type of boundary condition is often called a Neu-
mann boundary condition.. If
f approaches infinity and the derivative remains
bounded, then (2.3.6) implies
x
vxu
= x(O)= When the value of the function is
given at the boundary, this is often called the Dirichlet boundary condition.
© 2004 by Chapman & Hall/CRC
2.3. COOLING FIN AND TRIDIAGONAL MATRICES 71
2.3.3 Model
The above derivation is useful because (2.3.1) and (2.3.2) suggest a way to
discretize the continuous model. Let
x
l
be an approximation of x(lk) where
k = O@q. Approximate the derivative x
{
(lk + k@2) by (x
l+1
x
l
)@k. Then
equations (2.3.2) and (2.3.3) yield the finite di
erence approximation, a discrete
model, of the continuum model (2.3.4)-(2.3.6).
Let
x
0
be given and let 1 l ? q:
[N(x
l+1
x
l
)@k N(x
l
x
l1
)@k] + kFx
l
= ki(lk)= (2.3.7)
Let
l = q:
[f(x
vxu
x
q
) N(x
q
x
q1
)@k] + (k@2)Fx
q
= (k@2)i(qk)= (2.3.8)
The discrete system (2.3.7) and (2.3.8) may be written in matrix form.
For ease of notation we let
q = 4, multiply (2.3.7) by k and (2.3.8) by 2k,
E 2N + k
2
F so that there are 4 equations and 4 unknowns:
Ex
1
Nx
2
= k
2
i
1
+ Nx
0
>
Nx
1
+ Ex
2
Nx
3
= k
2
i
2
>
Nx
2
+ Ex
3
Nx
4
= k
2
i
3
and
2Nx
3
+ (E + 2kf)x
4
= k
2
i
4
+ 2fkx
vxu
=
The matrix form of this is DX = I where D is, in general, q × q matrix and X
and I are q × 1 column vectors. For q = 4 we have
D =
5
9
9
7
E N
0 0
N E N 0
0 N E N
0 0 2N E + 2fk
6
:
:
8
where X =
5
9
9
7
x
1
x
2
x
3
x
4
6
:
:
8
and I =
5
9
9
7
k
2
i
1
+ Nx
0
k
2
i
2
k
2
i
3
k
2
i
4
+ 2fkx
vxu
6
:
:
8
=
2.3.4 Method
The solution can be obtained by either using the tridiagonal (Thomas) algo-
rithm, or using a solver that is provided with your computer software. Let us
consider the tridiagonal system D{ = g where D is an q × q matrix and { and g
are q×1 column vectors. We assume the matrix D has components as indicated
in
D =
5
9
9
7
d
1
f
1
0 0
e
2
d
2
f
2
0
0
e
3
d
3
f
3
0 0 e
4
d
4
6
:
:
8
=
© 2004 by Chapman & Hall/CRC
72 CHAPTER 2. STEADY STATE DISCRETE MODELS
In previous sections we used the Gaussian elimination algorithm, and we
noted the matrix could be factored into two matrices
D = OX. Assume D is
tridiagonal so that
O has nonzero components only in its diagonal and subdi-
agonal, and X has nonzero components only in its diagonal and superdiagonal.
For the above 4 × 4 matrix this is
5
9
9
7
d
1
f
1
0 0
e
2
d
2
f
2
0
0
e
3
d
3
f
3
0 0 e
4
d
4
6
:
:
8
=
5
9
9
7
1
0 0 0
e
2
2
0 0
0
e
3
3
0
0 0 e
4
4
6
:
:
8
5
9
9
7
1
1
0 0
0 1
2
0
0 0 1
3
0 0 0 1
6
:
:
8
=
The p lan of action is (i) solve for
l
and
l
in terms of d
l
, e
l
and f
l
by
matching components in the above matrix equation, (ii) solve
O| = g and (iii)
solve X{ = |.
Step (i): For l = 1, d
1
=
1
and f
1
=
1
1
. So,
1
= d
1
and
1
= f
1
@d
1
.
For 2
l q 1, d
l
= e
l
l1
+
l
and f
l
=
l
l
= So,
l
= d
l
e
l
l1
and
l
= f
l
@
l
. For l = q, d
q
= e
q
q1
+
q
= So,
q
= d
q
e
q
q1
= These steps
can be executed provided the
l
are not zero or too close to zero!
Step (ii): Solve
O| = g.
|
1
= g
1
@
1
and for l = 2> ===> q |
l
= (g
l
e
l
|
l1
)@
l
.
Step (iii): Solve
X{ = |.
{
q
= |
q
and for l = q 1> ===> 1 {
l
= |
l
l
{
l+1
.
The loops for steps (i) and (ii) can be combined to form the following very
imp ortant algorithm.
Tridiagonal Algorithm.
(1) = a(1), (1) = c(1)/a(1) and y(1) = d(1)/a(1)
for i = 2, n
(i) = a(i)- b(i)*(i-1)
(i) = c(i)/(i)
y(i) = (d(i) - b(i)*y(i-1))/
(i)
endloop
x(n) = y(n)
for i = n - 1,1
x(i) = y(i) -
(i)*x(i+1)
endloop.
2.3.5 Implementation
In this section we use a MATLAB user defined function trid.m and the tridi-
agonal algorithm to solve the finite di
erence equations in (2.3.7) and (2.3.8).
The function
wulg(q> d> e> f> g) has input q and the column vectors d> e> f. The
output is the solution of the tridiagonal algebraic system. In the M
ATLA B code
fin1d.m lines 7-20 enter the basic data for the cooling fin. Lines 24-34 define
the column vectors in the variable list for trid.m. Line 38 is the call to trid.m.
© 2004 by Chapman & Hall/CRC
2.3. COOLING FIN AND TRIDIAGONAL MATRICES 73
the heat balance is computed in lines 46-54. Essentially, this checks to see if
the heat entering from the hot mass is equal to the heat lost o
the lateral
and tip areas of the fin. More detail about this will be given later. In the
trid.m function code lines 8-12 do the forward sweep where the OX factors are
computed and the O| = g solve is done. Lines 13-16 do the backward sweep to
solve X{ = |.
MATLAB Codes fin1d.m and trid.m
1. % This is a mo del for the steady state cooling fin.
2. % Assume heat di
uses in only one direction.
3. % The resulting algebraic system is solved by trid.m.
4. %
5. % Fin Data.
6. %
7. clear
8. n = 40
9. cond = .001;
10. csur = .001;
11. usur = 70.;
12. uleft = 160.;
13. T = .1;
14. W = 10.;
15. L = 1.;
16. h = L/n;
17. CC = csur*2.*(W+T)/(T*W);
18. for i = 1:n
19. x(i) = h*i;
20. end
21. %
22. % Define Tridiagonal Matrix
23. %
24. for i = 1:n-1
25. a(i) = 2*cond+h*h*CC;
26. b(i) = -cond;
27. c(i) = -cond;
28. d(i) = h*h*CC*usur;
29. end
30. d(1) = d(1) + cond*uleft;
31. a(n) = 2.*cond + h*h*CC + 2.*h*csur;
32. b(n) = -2.*cond;
33. d(n) = h*h*CC*usur + 2.*csur*usur*h;
34. c(n) = 0.0;
35. %
36. % Execute Tridiagonal Algorithm
© 2004 by Chapman & Hall/CRC
The output can be given as a table, see line 44, or as a graph, see line 55. Also,
74 CHAPTER 2. STEADY STATE DISCRETE MODELS
37. %
38. u = trid(n,a,b,c,d)
39. %
40. % Output as a Table or Graph
41. %
42. u = [uleft u];
43. x = [0 x];
44. % [x u];
45. % Heat entering left side of fin from hot mass
46. heatenter = T*W*cond*(u(2)-u(1))/h
47. heatouttip = T*W*csur*(usur-u(n+1));
48. heatoutlat =h*(2*T+2*W)*csur*(usur-u(1))/2;
49. for i=2:n
50. heatoutlat=heatoutlat+h*(2*T+2*W)*csur*(usur-u(i));
51. end
52. heatoutlat=heatoutlat+h*(2*T+2*W)*csur*(usur-u(n+1))/2;
53. heatout = heatouttip + heatoutlat
54. errorinheat = heatenter-heatout
55. plot(x,u)
1. function x = trid(n,a,b,c,d)
2. alpha = zeros(n,1);
3. gamma = zeros(n,1);
4. y = zeros(n,1);
5. alpha(1) = a(1);
6. gamma(1) = c(1)/alpha(1);
7. y(1) = d(1)/alpha(1);
8. for i = 2:n
9. alpha(i) = a(i) - b(i)*gamma(i-1);
10. gamma(i) = c(i)/alpha(i);
11. y(i) = (d(i) - b(i)*y(i-1))/alpha(i);
12. end
13. x(n) = y(n);
14. for i = n-1:-1:1
15. x(i) = y(i) - gamma(i)*x(i+1);
16. end
f = fvxu in (2.3.4) and (2.3.6). For larger f the solution or temperature should
be cl oser to the surrounding temperature, 70. Also, for larger
f the derivative
at the left boundary is very large, and this indicates, via the Fourier heat law,
that a large amount of heat is flowing from the hot mass into the right side of
the fin. The heat entering the fin from the left should equal the heat leaving
the fin through the lateral sides and the right tip; this is called heat balance.
© 2004 by Chapman & Hall/CRC
In Figure 2.3.2 the graphs of temperature versus space are given for variable
2.3. COOLING FIN AND TRIDIAGONAL MATRICES 75
Figure 2.3.2: Temperature for c = .1, .01, .001, .0001
2.3.6 Assessment
In the derivation of the model for the fin we made several assumptions. If the
thickness W of the fin is too large, there will be a varying temperature with the
vertical coordinate. By assuming the W parameter is large, one can neglect
any end e
ects on the temperature of the fin. Another problem arises if the
temperature varies over a large range in which case the thermal conductivity
N will be temp erature dependent. We will return to these problems.
Once the continuum model is agreed upon and the finite di
erence approxi-
mation is formed, one mu st be concerned about an appropriate mesh size. Here
an analysis much the same as in the previous chapter can be given. In more
complicated problems several computations with decreasing mesh sizes are done
until little variation in the numerical solutions is observed.
Another test for correctness of the mesh size and the model is to compute
the heat balance based on the computations. The heat balance simply s tates
the heat entering from the hot mass must equal the heat leaving through the
fin. One can derive a formula for this based on the steady state continuum
model (2.3.4)-(2.3.6). Integrate both sides of (2.3.4) to give
Z
O
0
0g{ =
Z
O
0
((2Z + 2W )@(W Z )f(x
vxu
x) + (Nx
{
)
{
)g{
0 =
Z
O
0
((2Z + 2W )@(W Z )f(x
vxu
x))g{ + Nx
{
(O) Nx
{
(0)=
Next use the boundary condition (2.3.6) and s olve for Nx
{
(0)
© 2004 by Chapman & Hall/CRC
76 CHAPTER 2. STEADY STATE DISCRETE MODELS
Nx
{
(0) =
Z
O
0
((2Z + 2W )@(W Z )f(x
vxu
x))g{
+f(x
vxu
x(O)) (2.3.9)
In the M
ATLA B code fin1d.m lines 46-54 approximate both sides of (2.3.9)
where the integration is done by the trapezoid rule and both sides are multiplied
by the cross section area, W Z . A large d ierence in these two calculations
indicates significant numerical errors. For q = 40 and smaller f = =0001, the
dierence was small and equaled 0=0023. For q = 40 and large f = =1, the
di
erence was about 50% of the approximate heat loss from the fin! However,
larger
q significantly reduces this dierence, for example when q = 320 and
large f = =1, then heat_enter = 3=7709, heat_out = 4=0550
The tridiagonal algorithm is not always applicable. Di
!culties will arise if
the
l
are zero or near zero. The following theorem gives conditions on the
components of the tridiagonal matrix so that the tridiagonal algorithm works
very well.
Theorem 2.3.1 (Existence and Stability) Consider the tridiagonal algebraic
system. If |
d
1
| A |f
1
| A 0, |d
l
| A |e
l
| + |f
l
|, f
l
6= 0, e
l
6= 0 and 1 ? l ? q,
|d
q
| A |f
q
| A 0> then
1. 0
? |d
l
| |e
l
| ? |
l
| ? |d
l
| + |e
l
| for 1 l q (avoids division by small
numbers) and
2. |
l
| ? 1 for 1 l q (the stability in the backward solve loop).
Proof. The proof uses mathematical induction on
q. Set l = 1: e
1
= 0 and
|
1
| = |d
1
| A 0 and |
1
| = |f
1
|@|d
1
| ? 1=
Set i A 1 and assume it is true for l 1:
l
= d
l
e
l
l1
and
l
= f
l
@
l
.
So,
d
l
= e
l
l1
+
l
and |d
l
| |e
l
||
l1
| + |
l
| ? |e
l
|1 + |
l
|= Then |
l
| A
|d
l
| |e
l
| |f
l
| A 0= Also, |
l
| = |d
l
e
l
l1
| |d
l
| + |e
l
||
l1
| ? |d
l
| + |e
l
|1=
|
l
| = |f
l
|@|
l
| ? |f
l
|@(|d
l
| |e
l
|) 1=
2.3.7 Exercises
1. By hand do the tridiagonal algorithm for 3{
1
{
2
= 1> {
1
+4{
2
{
3
= 2
and
{
2
+ 2{
3
= 3.
2. Show that the tridiagonal algorithm fails for the following problem
{
1
{
2
= 1> {
1
+ 2{
2
{
3
= 1 and {
2
+ {
3
= 1.
3. In the derivation of the tridiagonal algorithm we combined some of the
loops. Justify this.
4.
ment with di
erent values of W = =05> =10> =15 and .20. Explain your results and
evaluate the accuracy of the model.
© 2004 by Chapman & Hall/CRC
Use the code fin1d.m and verify the calculationsinFigure 2.3.2. Experi-
2.4. SCHUR COMPLEMENT 77
5. Find the exact solution of the fin problem and experiment with di
erent
mesh sizes by using
q = 10> 20> 40 and 80. Observe convergence of the discrete
solution to the continuum solution. Examine the heat balance calculations.
6. Modify the above model and code for a tapered fin where
W = =2(1 {)+
=1{.
7. Consider the s teady state axially symmetric heat conduction problem
0 =
ui + (Nux
u
)
u
, x(u
0
) = jlyhq and x(U
0
) = jlyhq. Assume 0 ? u
0
? U
0
.
Find a discrete model and the solution to the resulting algebraic problems.
2.4 Schur Complement
2.4.1 Intro duction
In this section we will continue to discuss Gaussian elimination for the solution
of D{ = g. Here we will examine a block version of Gaussian elimination. This
is particularly useful for two reasons. First, this allows for e
!cient use of the
computer’s memory hierarchy. Second, when the algebraic equation evolves
from models of physical objects, then the decomposition of the object may
match with the blocks in the matrix
D. We will illustrate this for steady state
heat diusion models with one and two space variables, and later f or models
with three space variables.
2.4.2 Applied Area
In the previous section we discussed the steady state model of diusion of heat
in a cooling fin. The continuous model has the form of an ordinary di
erential
equation with given temperature at the boundary that joins the hot mass. If
there is heat diusion in two directions, then the model will be more compli-
cated, which will be more carefully described in the next chapter. The objective
is to solve the resulting algebraic system of equations for the approximate tem-
perature as a function of more than one space variable.
2.4.3 Model
The continuous models for steady state heat diusion are a consequence of the
Fourier heat law applied to the directions of heat flow. For simplicity assume the
temperature is given on all parts of the boundary. More details are presented
in Chapter 4.2 where the steady state cooling fin model for diusion in two
directions is derived.
Continuous Models:
Di
usion in 1D. Let x = x({) = temperature on an interval.
0 =
i + (Nx
{
)
{
and (2.4.1)
x(0)> x(O) = given= (2.4.2)
© 2004 by Chapman & Hall/CRC