T h e right-hand side becomes
lim/1„ [ l _ _ * * « . ] . (2.6.28) I f hx and ht are related such that ht/hx2 > ^, then the argument in
(2.6.28) is of magnitude greater than unity. T h e logarithm is positive (plus a phase factor) and hence the solution will diverge, thereby entirely failing to represent the desired solution. I n order to guarantee the con
vergence of the approximation, w e must require
\ < \ (2-6.29) even as ht,hx-+ 0.
Requirements such as (2.6.29) are very frequent in numerical studies, and w e shall encounter such relations often in our later analysis.
2.7 Matrix Form of Difference Equations
In our later work w e shall frequently use a compact matrix notation for the simultaneous equations resulting from difference approximation to differential equations. T o illustrate the construction of the matrix form of the equation, w e n o w consider the t w o examples of partial difference equations given in Section 2.5.
W e approximate the heat-flow equation by the simple difference relation
Tk,i+i — Tktl — -j— [Tk+ltl — 2Tktl + Tk_ia], (2.7.1)
X
W e assume the boundary conditions as before, Eqs. (2.5.3a,b). Equation (2.7.1) can be factored in the form
h
T h e above set of simultaneous equations can evidently be written
- α 1 0 ... 0 0"
1 - α 1 ... 0 0
0 0 0 ... -OL 1 0 0 0 ... 1 -oc
L^V-I.ZJ
L^tf-l.l+lJ
ι,ι+ι
^2,1+1
(2.7.4) I f w e denote the matrix in (2.7.4) as A , and define the vector ψ , as
T2l
(2.7.5)
then the entire set of equations can be written
hi Α ψ( = ψί + 1. (2.7.6)
I f w e denote the starting vector as ψ0 , then E q . (2.7.6) can also be written
[(Α,/ΑΪΪΑΐ'ψο = ψ , .
N o t e that the matrix A is of tridiagonal form.
For the Laplace difference equation, w e have the relation
h2 ^ h2
(2.7.7)
(2.7.8)
By simple factoring, w e have
r2tyk+lJ +Φΐ€-υ) — βΦν + ΦΜ+1 + ΦkJ-l = 0,
where r2 - h2\h,2 and β = 2(1 + r2).
(2.7.9)
2.7 M A T R I X F O R M O F D I F F E R E N C E E Q U A T I O N S
W e now define the extended vector ψ
Ψ =
T h e set of equations (2.7.11) can then be written Β ψ = 0 ,
Each of the elements of Β is a (Κ — 1) by (Κ — 1) square submatrix.
A matrix of the form of Β is a block tridiagonal matrix. T h e procedure can obviously extended to difference equations with more unknowns and also to equations with variable coefficients and/or variable mesh spacing.
References
There are several books devoted to the calculus of finite differences. Parti
cularly exhaustive treatments are included in references 7, 2, and 3. Somewhat shorter treatments are found in 4 and 5. Reference 5 also includes many other aspects of the difference calculus. T h e application of finite differences is not limited to approximating differential equations. Other uses for the calculus are described very cogently in 4 and 6.
L Jordan, C , "Calculus of Finite Differences." Chelsea, N e w York, 1947.
2. Milne-Thompson, L. M . , " T h e Calculus of Finite Differences.'' Macmillan, London, 1933.
3. Fort, T., "Finite Differences and Difference Equations in the Real Domain.'' Oxford Univ. Press, London and N e w York, 1948.
4. Hildebrand, F. B., "Introduction to Numerical Analysis." McGraw-Hill, N e w York, 1956.
5. Hildebrand, F. B., "Methods of Applied Mathematics." Prentice-Hall, Englewood Cliffs, N e w Jersey, 1952.
6. Hamming, R. W . , "Numerical Methods for Scientists and Engineers." M c G r a w
4. Prove that the differentiation operator d/dx may be formally written
— = τ sinh 1
-dx h
2
and then develop expansion (2.3.17).
P R O B L E M S 8 3
5. Find the solution to the following:
(a) A*gk + Agk - gk = 0, (b) g(k + 3) - 8*(Λ) = 0.
6. Find the solutions to the eigenvalue problem + A«y* = 0, 0
g(Q) = l,
with boundary conditions
4Vo = 0; Vy* = 0.
Show the convergence to the analytic solution.
A n infinite medium, void of neutrons, has a constant scattering cross-section and a constant fission cross-section. There is no capture. T h e fission process produces 1 prompt neutron and another neutron delayed exactly τ sees.
T h e total diffusion time of any neutron is exactly τ sees. One neutron is introduced at time t = 0 into the assembly.
[a] Derive an expression for the neutron population as a function of time.
[b] T h e assembly is to be scrammed when more than 1015 neutrons are present. W h e n should scram occur ?
8. Approximate the wave equation 3*φ dt2 with the difference approximation
θ2φ
82φ
Solve the difference equation and compare the solutions with the analytic solution to the differential equation. Under what circumstances will the difference equation converge ?
9. The determinant, ΔΝ, of the Ν by Ν matrix - « 1 0 ·
1 - a 1 0
arises from applying the operator δ2. By assuming A0 = 1, show the recurrence relation
is valid.
AN = - AN_2
10. Derive a difference approximation to the Laplacian operator in r, ζ coordi
nates which is accurate to order h\ and h\.
11. By noting that the series Σ^ο* einkn/K is geometric, prove
(a)
(b)
[The other orthogonality relations given in Section 2.5 may be developed in a similar manner].
12. The purpose of this problem is to illustrate the construction of a difference equation by physical rather than mathematical means. Heat is transferred from a cylindrical fuel element of radius R() encased in cladding of uniform thickness a, centered in a cylindrical coolant channel of radius RK . The coolant channel is assumed insulated at the outer boundary; the coolant mass flow rate and inlet temperature are assumed fixed. Calculate the tem
perature of the fuel, cladding, and coolant as functions of r and z> the radial and axial coordinates. Neglect any heating in the coolant or cladding as a result of neutron or photon irradiation. Assume thermal conductivities are independent of temperature in all materials; neglect axial heat flow; neglect any thermal resistance between cladding and fuel. Assume the heat source in the fuel depends only on axial position and that the coolant is completely turbulent.
(a) Write down the usual energy conservation equations and boundary conditions needed to solve the problem. Let η be the heat transfer coefficient between the cladding and the coolant, and kf the fuel thermal conductivity.
(b) Solve the resulting equations analytically for the temperature as a function of r and ζ in the fuel, cladding, and coolant.
(c) Divide the assembly into a number of annuli such that the radius tj of the 7th one is given by r, = jh. By performing energy balances directly on zones j — 1, j\ and j + 1, show that
-^i-
(Γ<_ι - T ) + ^) - ( Τί +1 - T,) =S (zW) - rU,
where T} is temperature at rt and S(z) is the heat source density. What happens at j = 1, i.e., how is the temperature of the centerline to be deter
mined ? Show that the above result is correct to first order.
(d) By using the cross-sectional area of flow at the average radius of an interval, show that
- ^ - ( T , - ! - T^r^ + r,) + - £ - ( Τ ,+1 - Τ, ) (Γ| + = S(z)(r) - rJ.J.
Eliminate r^x, rj, and ri+1 in terms of j and h. Show that the truncation error is of order h.
(e) As the next improvement, for the source use the volume between tj — h/2 and r} + h/2f and show that the heat flow equation is
- £ - [ ( T f- i - T^r^ + r,) + ( T ,+l - T,)(r, + r ,+ 1) ] = S ( * ) ( r J+i - rj.j).
P R O B L E M S 85