A Numerical Solution Technique for Linear
Boundary Value Problems of Ordinary Differential Equations
Abstract:
A method for the
solution of ordinary linear differential equations (OLDEs) is introduced. The
significance of this method is in its directness and simplicity. An Nth order
OLDE can be considered as a directional derivative functional made of a scalar
product of two vectors in an (N+1) dimensional space. One of these vectors
consists of (N+1) independent functions among the total (K+1) (K: number of
numerical integration steps) functions and the other is the gradient of the
solution at that step with respect to the components of this vector. With the
help of the recursion relation that is obtained by differencing the OLDE, the
gradient vectors can be obtained recursively, and numerically stored in digital
computers that are programmable through any of the present languages. In other
words, through this method we can effectively self start the solution as would
be in an initial value problem.
Introduction:
The differential
equations arise in scientific and engineering problems naturally as an
expression of the change of the independent variable with respect to the
control variable or variables. In the case of more than one control variable
the problem is called a partial differential equation. The full solution,
however is obtained only after the introduction of some conditions that the
desired solution or its derivatives have to satisfy at the ends. If these end
conditions are all at the one and same side, the problem is a self starter and
properly called an initial value problem. For these types of problems both
numerical and, if available, analytical solutions are found straightaway. The
analytical solution, if available, is of course the first choice. However if
the solution cannot be made expressible in terms of any known analytical
functions, the only resort that remains is the numerical approach. For the more frequent type where the end
(boundary) conditions are given at the different ends, the self start is not
possible and the numerical solution is traditionally found from either, a
(1)shooting method or, a (2) matrix inversion method [1].
In this paper we
introduce a straight forward numerical solution technique for any order of OLDEs.
The linearity means that, each term of the equation can only contain the
dependent variable and its derivatives of any order to the first power. It
should be mentioned that, when the type of the problem permits a separation of
variables approach, some partial differential equations can be reduced to OLDEs
and hence they too can yield to our method directly, otherwise they can be
dealt with similarly as will be shown in a subsequent paper.
The well known
boundary value problems in physical and engineering applications are Laplace's
(or similarly Poisson's: for example, Uxx+Uyy+Uzz=C(x,y,z,t),
for the unknown function U(x,y,z,t), where subscripts denote partial
differentiation with respect to the denoted variable as many times as the
repetition, an elliptic partial differential equation (PDE)) equation, wave
equation (for example, Uxx+Utt=0, a hyperbolic PDE),
diffusion equation (or heat equation, for example, Uxx+Uyy+Uzz==Ut,
a parabolic PDE) [2], and the telegram equation. Additionally, there is an
integro-differential variety, known as radiative transfer equation, which under
suitable approximations can be converted to a diffusion type equation [3].
For clarity of
demonstration and to be more specific, we will show only the application of our
technique to two-stream approximation solution of the radiative transfer
equation which is a second order OLDE with two boundary conditions (each at
different ends, of course). Again it should be stressed that the method is
general and applies to any order of differential equations whether they are
ordinary or partial as long as linearity is satisfied. From this point on,
whenever a reference is made to a differential equation, it is meant to be an
OLDE type with the prescribed boundary conditions.
There is another
variation of the two-stream approximation that is known as delta-Eddington
method [3; 4]. In either case the final forms of the governing differential
equations are:
Y''(x) + f(x).Y'(x) +
g(x).Y(x) = h(x) ....................................................................................(1)
and
Z''(x) + j(x).Z'(x) +
k(x).Z(x) = m(x)
.....................................................................................(2)
where primes denote
differentiation with respect to x, and f(x), g(x), h(x), j(x), k(x), m(x) are
the variable coefficients for OLDEs of the sought after solutions Y(x) and
Z(x). They result from de-coupling two linear first order differential
equations. Alternatively these equations can be cast into a form without the
first order differentials through a proper change of the variables [5]. In the
following we shall demonstrate our method for both forms. If the functions f(x)
etc. are constants, the solutions can be found analytically as positive and
negative exponentials of x times some constants. The final solution can be
obtained after the application of the given boundary conditions. But if the
path the photons travel inside the medium has position dependent optical
properties (absorption and scattering coefficients times material densities),
then the functions f(x) etc. should vary with respect to x, hence they are
named solutions for inhomogeneous paths. To give an idea on the current state
of the inhomogeneous path modeling in the atmospheric applications, we only
need to mention that even though a measurement of the vertical inhomogeneity in
the clouds was reported by Roewe and Liou in 1978 [6], since that time the
effects of inhomogeneity within an emitting-scattering and absorbing medium,
even along the vertical direction has not been a wide area of study. In the
current state of affairs mostly the homogeneous path solutions are contended
with. The main reason for that we suspect, is the mathematical complexity of
the required solutions and the ensuing complications that would arise in the
remote sensing inversion algorithms in comparison to the present ones.
Ours and the previous methods:
In order to construct
the numerical solutions to the above equations (1) or (2), the domain x can be
divided into n layers of equal thickness d, through m=n+1 levels with labels
1,2,...n, m; this is differencing. For
example the first of the above differential equations in difference form gives,
[Y(i+2)-2.Y(i+1)+Y(i)]/(d2)
+ f(i+1).[Y(i+2)-Y(i)]/(2.d) + g(i+1).Y(i+1) =
h(i+1) ...................(3)
or
Y(i+2) = { [2-g(i+1).d2].Y(i+1)
+ [f(i+1).d/2-1].Y(i) + h(i+1).d2 } / [1+f(i+1).d/2].....................(4)
where array numbers
within the parentheses indicate the value of Y at that level. Each layer
extends to infinity horizontally and is horizontally homogenous, hence the name
plane parallel. If the end conditions were given at the same boundary, as say
Y(1) and Y'(1) (the derivative of Y at 1), then the solution would be obtained
through Euler’s integration scheme(or any more sophisticated one such as a
Runge-Kutta’s) as would be in an initial value problem. In two dimensions, when
only Y(1) and Y(m) are given, we have a true boundary value problem and without
the apriori knowledge of Y(2), literature assures us that we cannot propagate
the solution to the differential equation up to Y(m). Yet in this paper we will
demonstrate that we easily can.
If we were using a
shooting method, the numerical solution at one boundary would start with a
trial and error parameter Y(2) which, if correctly chosen would give the
boundary match at the other end of the boundary x(m). Sometimes in practice
Y(m) that is obtained from the shooting method might wildly swing around the
boundary value at x(m) even for small variations in the trial value Y(2). In
our opinion, the best and really required place to use a trial and error method
would be when the differential equation is nonlinear to begin with.
On the other hand if
we were using a matrix inversion method, the differential equation would be
solved for n thin layers of the medium in which the optical properties and so
the functions f(x) etc. would be taken to be constants. The resultant set of n
equations are matched at their common boundaries and (2n-2) redundant
coefficients are needed to be found before the desired solution is obtained.
This corresponds to the traditional numerical solution of partial differential
boundary value problems where the unknown are all the nodal values.
In atmospheric
radiative transfer computations in a plane parallel medium, a package named
DISORT (DIScrete Ordinate Radiative Transfer) [7] solves an n-set of second
order OLDEs with s streams by matrix inversion. The size of the matrix is
(s.n)-by-(s.n), where s is the number of streams (two for a two stream solution
[4;7]), whereas in our method it is only s-by-s.
Thus one of the
advantages of our method is relief from large matrix inversions and the other
advantage is the speed gained through the directness of the method because it
is free of redundant calculations for the intermediate boundaries. It appears
that these existing methods in the literature have not yet taken into account
the full capabilities of the digital computers.
If we return to the
recursive relation (3), where s=2 and the number of layers is n, we see that we
have m=n+1 levels corresponding Y(1),Y(2), ..,Y(m). Also we observe that any of
these Y(i)'s depends on Y(1), Y(2) and H(2)=h(2).d2/[1+f(2).d/2] for i greater than 2. Y(1), Y(2) and H(2)
are independent of each other and so we can choose three orthonormal vectors as
a basis for a three dimensional vector space as,
e1 = ÑY(1) = (1,0,0), e2 = ÑY(2) = (1,0,0), and e3 = ÑP(2) = (1,0,0).
Also note that H(i)=
H(i-1) . [H(i)/H(i-1)]. The ratio
[H(i)/H(i-1)] =
H(i)/H(i-1). H(i-1)/H(i-2). H(i-2)/H(i-3)...... H(3)/H(2) is independent
of H(i) or H(2), since H(3)/H(2) =
[H(2) + DH(2)]
/ H(2) = 1 + D
ln[H(2)] and so forth.
Then we can rewrite
the relation (3) which is for i=1, as a directional derivative of Y(3) along
the vector [Y(1) e1 + Y(2) e2
+ H(2) e3],
Y(3) = [Y(1) e1+Y(2) e2+ H(2) e3].ÑY(3) or Y(3) = [Y(1), Y(2), H(2)].ÑY(3)
............(5)
where
ÑY(3)
= { -[f(2)/(2.d)-1/(d.d)],
[2/(d.d)-g(2)], 1 } /
[1/(d.d)+f(2)] ........................................(6)
Similar gradient
operations for all other Y(i)'s up to Y(m) are to be done by using the
recursive relations similar to (4) and (5). A gradient operation for Y(i)'s and
h(i-1)'s in left hand side vector of
(3) recursively until the gradient with respect to the basis of e1, e2 and e3 is
obtained. Along the way all the numerics get taken care of in the ÑY(i)
vector (or array for the digital computer) as in
Y(i) = [Y(1), Y(2),
h(2)].Ñ
[Y(i)].................................................................................................(7)
At the end of the
calculations we need to only to match Y(m) to the other end boundary value and
determine Y(2). With Y(2) found, the next step is to perform the scalar product
with each i step of the gradient vector and so find the total solution for all
the steps. In effect, the whole solution is obtained as if it were an initial value problem.
In the following we
will demonstrate the method with three examples. For this purpose, we shall
only deal with one of the above equations and its alternative forms as in:
Y''(x) - Y(x) = 0
......................................................................................................................(8)
Y''(x) + g(x).Y(x) =
h(x) ..........................................................................................................(9)
Y''(x) + f(x).Y'(x) +
g(x).Y(x) = h(x)
.......................................................................................(10)
The first one above
corresponds to the case where g(x) = -1, f(x) = h(x) = 0. This is the simplest case and the best one
that elucidates the simplicity of the method.
The second one above
has no Y(x) term (f(x) = 0), but otherwise equivalent to the third one. The
third one is the most general form whose method of attack has been just
mentioned above. It is a kind that can arise in the study of many other
problems in physics and engineering.
The example solutions for the three cases and discussions:
The solutions to the
sample function are given in the Fortran codes with their outputs. The method
could be coded in any programming language that allows matrix definitions and
algebraic operations. These languages include but not limited to Basic, C,
Pascal, and other mathematical packages such as Mathematica, Mathlab, IDL, etc.
The heart of the method lies in the vectorization of Y by gradient operation.
This technique bypasses the impossibility of doing mixed symbolic and numerical
computations at the same time with the above mentioned languages. Through the
vectors, the symbolics are carried through the mth value where the
upper boundary condition is satisfied and Y(2) is then found at once. The three
examples are given below with sufficient comments in them in order to clarify
the method. For higher precision, the resolution can be made finer whenever
required. The examples are for the test functions mentioned in the codes. Since
the sole purpose of this report is to present the method, no further work on
the convergence of the integration scheme or no research for the fastest
integration scheme is carried out here.
The examples:
1) The
simplest example that demonstrates our method, especially the vectorization
scheme is the following:
We used Y(x) =
exp(x)+exp(-x) as our test function. With only eleven steps between x1=0 and
xm=1, the percent errors are found to be less than 0.01\%. The differential
equation is
Y''(x) - Y(x) = 0,
where g(x) = -1. The equation is rewritten in difference form as
Y3-2Y2+Y1-Y1*d2=0,
where d = dx or the step size.
Y3 =
Y1*(1-d2)+Y2, Y4 = Y2*(1-d2)+Y3 =
Y2*(1-d2)+ Y1*(1-d)+Y2 = Y2*(2-d2)+ Y1*(1-d2) etc. Or in a more generel form
Yi=Yi-1*(1-d2)+Yi-2 =
Y1*a+Y2*b.
If we were given the
numerical value of Y2, the solution would be a self starter and we would obtain
all Yi's at once. For the starting values, both Y1 and Y2 are vectorized.
Y(1,1) = ÑY(1).e1, Y(1,2) = ÑY(1).e2, Y(2,1) = ÑY(2).e1, Y(2,2) = ÑY(1).e2
and these
initialization values as defined in the data statement. d2 is (dx)2, where dx = (1-0)/10, since
we divided the integration interval into ten equal steps.
a =Y(0) and b =Y(1)
are the boundary values of the function Y. A close inspection of the Fortran
code below reveals the method and shows its surprisingly economical logics.
Code 4 further reduces
the dimensionality of the solution by absorbing the inhomogeneous part into the
second vector component.
.........................................................................................................................................................
Code 1:
An example for a
linear homogeneous ordinary differential equation with constant coefficients,
and with no Y(x) term:
c Example: Y''(x) - Y(x) = 0
real g(11), Y(11,2), Yx(11) !
n=10, m=n+1=11
data/Y(1,1),Y(1,2),Y(2,1),Y(2,2),g,d2,a/1,0,0,1,11*-1,.01,2/
b = exp(1)+exp(-1) ! the upper
boundary
do i=2, 10, 1
Y(i+1,1) = (2-d2*g(i))*Y(i,1) -
Y(i-1,1)
Y(i+1,2) = (2-d2*g(i))*Y(i,2)
-Y(i-1,2)
end do
c=b-Y(m,1)*a/Y(m,2) from
b=Ym=Y(m,1)*a+Y(m,2)*c
c: print*,' x Ysoln
Ytest % error'
do i=1,11,1
Yx(i) =
Y(i,1)*a+Y(i,2)*(b-Y(11,1)*a)/Y(11,2)
Ytest =
exp((i-1)/10.)+exp(-(i-1)/10.)
print*,(i-1)/10.,Yx(i),Ytest,(Yx(i)/Ytest-1)*100
end do
end
The output:
x Ysoln
Ytest % error
........................................................................................
0.000000E+00 2.000000
2.000000 0.000000E+00
1.000000E-01 2.010083 2.010008 3.736395E-03
2.000000E-01 2.040267 2.040133 6.567768E-03
3.000000E-01 2.090855 2.090677 8.495901E-03
4.000000E-01 2.162350 2.162145 9.494203E-03
5.000000E-01 2.255469 2.255252 9.620252E-03
6.000000E-01 2.371142 2.370930 8.939702E-03
7.000000E-01 2.510527 2.510338 7.531493E-03
8.000000E-01 2.675018 2.674870 5.535146E-03
9.000000E-01 2.866258 2.866173 2.961336E-03
1.000000E+00 3.086161
3.086161 0.000000E+00
.........................................................................................................................................................
Code 2:
c An example for
non-homogeneous ordinary linear differential equation
c with variable
coefficients, with no Y(x) term (f(x)=0), Y" + g(x)*Y = h(x)
parameter (m=11) ! m = number of
solution intervals + 1
real x(m), Y(m,3), Yx(m), f(m),
g(m)
x1 = 0 ! The starting point of the solution.
xm = 1 !The end point of the
solution.
Y1 = exp(x1)+exp( x1) ! Boundary condition 1 for Y at x=x1.
Ym = exp(xm)+exp(-xm) ! Boundary
condition 2 for Y at x=xm.
do i = 1, m ! Obtain x, g(x), and
h(x).
x(i) = x1 + (i-1.)*(xm-x1) /
(m-1.)
g(i) = x(i)
h(i) = (1+x(i))*(exp(x(i)) +
exp(-x(i)))
end do
d = (xm-x1) / (m-1.) ! Solution step
size
c d2 = d*d ! In case needed.
Y(1,1)=1
Y(1,2)=0
Y(1,3)=0
Y(2,1)=0
Y(2,2)=1
Y(2,3)=0
do i = 2, m-1 ! Iterative solution
as vector components.
Y(i+1,1) = (2-d2*h(i))*Y(i,1) -
Y(i-1,1)
Y(i+1,2) = (2-d2*h(i))*Y(i,2) -
Y(i-1,2)
Y(i+1,3) = (2-d2*h(i))*Y(i,3) -
Y(i-1,3) + h(i)*d2
end do
c BC2:Ym = Y(m,1)*Y1 +
Y(m,2)*Y2 + (Ypm=Y(n,3)), thus the unknown Y2 is
Y2 = (Ym-Y(m,1)*Y1-Y(m,3)) / Y(m,2)
print*,' x
Ytest Ysoln % error'
do i = 1, m ! Construct Ysoln as
Yx(i) and compare with Ytest.
Yx(i) = Y(i,1)*Y1 + Y(i,2)*Y2 +
Y(i,3)
Ytest = exp(x(i)) + exp(-x(i))
perror = (Ytest-Yx(i)) / Ytest *
100. ! Percent error.
print*, x(i), Ytest, Yx(i), perror
end do
stop
end
The output:
x Ytest
Ysoln % error
..........................................................................................
0.000000E+00 2.000000 2.000000
0.000000E+00
1.000000E-01 2.010008 2.010094
-4.270101E-03
2.000000E-01 2.040133 2.040288 -7.582713E-03
3.000000E-01 2.090677 2.090883
-9.863766E-03
4.000000E-01 2.162145 2.162384
-1.107803E-02
5.000000E-01 2.255252 2.255506
-1.124626E-02
6.000000E-01 2.370930 2.371178
-1.045673E-02
7.000000E-01 2.510338 2.510559
-8.797637E-03
8.000000E-01 2.674870 2.675042
-6.421313E-03
9.000000E-01 2.866173 2.866271
-3.437876E-03
1.000000E+00 3.086161 3.086162
-1.886656E-05
.........................................................................................................................................................
Code 3:
c An example for a
second order non-homogeneous ordinary
c linear differential
equation with variable coefficients:
c Y" + f(x)*Y' +
g(x)*Y = h(x)
parameter (m=11) ! m = number of
solution intervals + 1
real Y(m,3), Yx(m), g, gx(m), h,
hx(m), f, fx(m)
Ytest(x) = exp(x) + exp(-x)
f(x) = x
g(x) = x
h(x) = Ytest(x)*(1+g(x)) +
f(x)*(exp(x)-exp(-x))
x1 = 0 ! The start point of the solution.
xm = 1 ! The end point of the solution.
Y1 = Ytest(x1) ! Boundary condition 1 (BC1) for Y at x=x1.
Ym = Ytest(xm) ! Boundary
condition 2 (BC2) for Y at x=xm.
do i = 1, m ! Obtain x, f(x),
g(x), and h(x) at i.
x = x1 + (i-1.)*(xm-x1) / (m-1.)
fx(i) = f(x)
gx(i) = g(x)
hx(i) = h(x)
end do
d = (xm-x1) / (m-1.) ! Solution
step size
c d2 = d*d ! In case needed.
Y(1,1) = 1
Y(1,2) = 0
Y(1,3) = 0
Y(2,1) = 0
Y(2,2) = 1
Y(2,3) = 0
do i = 2, m-1 ! Iterative solution
as vector components.
Y(i+1,1)=(2-d2*gx(i))/(1.+fx(i)*d/2.)* Y(i,1) +
&
(fx(i)*d/2.-1.)/(1.+fx(i)*d/2.)*Y(i-1,1)
Y(i+1,2)=(2-d2*gx(i))/(1.+fx(i)*d/2.)* Y(i,2) +
& (fx(i)*d/2.-1.)/(1.+fx(i)*d/2.)*Y(i-1,2)
Y(i+1,3)=(2-d2*gx(i))/(1.+fx(i)*d/2.)* Y(i,3) +
&
(fx(i)*d/2.-1.)/(1.+fx(i)*d/2.)*Y(i-1,3) +
& (hx(i)*d2) /(1.+fx(i)*d/2.)
end do
c BC2:Ym = Y(m,1)*Y1 +
Y(m,2)*Y2 + (Ypm=Y(m,3)), thus the unknown Y2 is
Y2 = (Ym-Y(m,1)*Y1-Y(m,3)) /
Y(m,2)
print*,' x
Ytest Ysoln % error'
do i = 1, n ! Construct Ysoln as
Yx(i) and compare with Ytest.
x = x1 + (i-1.)*(xm-x1) / (m-1.)
Yx(i) = Y(i,1)*Y1 + Y(i,2)*Y2 +
Y(i,3)
perror = (Ytest(x)-Yx(i)) /
Ytest(x) *100. ! Percent error.
print*, x, Ytest(x), Yx(i),
perror
end do
stop
end
The output:
x
Ytest Ysoln % error
..........................................................................................
0.000000E+00 2.000000 2.000000 0.000000E+00
1.000000E-01 2.010008 2.010127 -5.883275E-03
2.000000E-01 2.040133 2.040352
-1.069130E-02
3.000000E-01 2.090677 2.090974
-1.420865E-02
4.000000E-01 2.162145 2.162498
-1.631583E-02
5.000000E-01 2.255252 2.255634
-1.692326E-02
6.000000E-01 2.370930 2.371310
-1.602770E-02
7.000000E-01 2.510338 2.510682
-1.371733E-02
8.000000E-01 2.674870 2.675142
-1.015598E-02
9.000000E-01 2.866173 2.866331
-5.509147E-03
1.000000E+00 3.086161 3.086162
-2.659197E-05
.........................................................................................................................................................
Code 4:
c An example for a second order non-homogeneous ordinary
c linear differential equation with variable coefficients:
c Y" + f(x)*Y' + g(x)*Y = h(x)
c234567
parameter (m=11) ! m =
number of solution intervals + 1
real Y(m,2), Yx(m), g,
gx(m), h, hx(m), f, fx(m)
Ytest(x) = exp(x) + exp(-x)
f(x) = x
g(x) = x
h(x) = Ytest(x)*(1+g(x)) + f(x)*(exp(x)-exp(-x))
x1 = 0 ! The start point of the solution.
xm = 1 ! The end point of the solution.
Y1 = Ytest(x1) ! Boundary condition 1 (BC1) for Y at x=x1.
Ym = Ytest(xm) ! Boundary condition
2 (BC2) for Y at x=xm.
do i = 1, m ! Obtain x,
f(x), g(x), and h(x) at i.
x = x1 + (i-1.)*(xm-x1) /
(m-1.)
fx(i) =
f(x)
gx(i) =
g(x)
hx(i) =
h(x)
end do
d = (xm-x1) / (m-1.) !
Solution step size
d2 = d*d ! In case needed.
Y(1,1) = Y1
Y(1,2) = 0
Y(2,1) = 0
Y(2,2) = 1
c do i = 2, m-1
c Y(i+1,1) =
Y(i-1,1)*(fx(i)*d/2.-1.)+(hx(i)*d2)
c end do
do i = 2, m-1 ! Iterative
solution as vector components.
Y(i+1,1)=(2-d2*gx(i))/(1.+fx(i)*d/2.)* Y(i,1) +
&
(fx(i)*d/2.-1.)/(1.+fx(i)*d/2.)*Y(i-1,1) +
& (hx(i)*d2)
/(1.+fx(i)*d/2.)
Y(i+1,2)=(2-d2*gx(i))/(1.+fx(i)*d/2.)* Y(i,2) +
&
(fx(i)*d/2.-1.)/(1.+fx(i)*d/2.)*Y(i-1,2)
end do
c BC2:Ym = Y(m,1) + Y(m,2)*Y2 + (Ypm=Y(m,3)), thus the unknown Y2 is
Y2 = (Ym-Y(m,1)) / Y(m,2)
print*,' x
Ytest Ysoln % error'
do i = 1, m ! Construct
Ysoln as Yx(i) and compare with Ytest.
x = x1 + (i-1.)*(xm-x1) /
(m-1.)
Yx(i) = Y(i,1) + Y(i,2)*Y2
perror = (Ytest(x)-Yx(i)) /
Ytest(x) *100. ! Percent error.
print*, x, Ytest(x), Yx(i),
perror
end do
stop
end
c----------The Output-------------------------------------------------
c x Ytest Ysoln % error
c 0.000000E+00 2.000000 2.000000
0.000000E+00
c 1.000000E-01 2.010008 2.010126
-5.856651E-03
c 2.000000E-01 2.040133 2.040351
-1.064341E-02
c 3.000000E-01 2.090677 2.090973
-1.414938E-02
c 4.000000E-01 2.162145 2.162496
-1.624821E-02
c 5.000000E-01 2.255252 2.255632
-1.685130E-02
c 6.000000E-01 2.370930 2.371309
-1.595621E-02
c 7.000000E-01 2.510338 2.510680
-1.364074E-02
c 8.000000E-01 2.674870 2.675140
-1.008760E-02
c 9.000000E-01 2.866173 2.866329
-5.450605E-03
c 1.000000 3.086161 3.086161
1.738217E-05
cStop - Program terminated.
References:
[1]- Hildebrand, F.
B., "Introduction to Numerical Analysis", 2nd Ed., 669 pages, Dover
Publications Inc., New York, USA, 1974.
[2]- James, M. L.,
Smith, G. M., Wolford, J. C., "Applied Numerical Methods for Digital
Computation", 4th Ed., 709 pages, HarperCollins College
Publishers, New York, USA, 1993.
[3]- Liou, Kuo-Nan,
"An Introduction to Atmospheric Radiation", 392 pages, Academic Press
Inc., San Diego, USA, 1980.
[4]- Lenoble, J.,
"Radiative Transfer in Scattering and Absorbing Atmospheres: Standard
Computational Procedures", 300 pages, Virginia USA, A. Deepak
Publishing, 1985.
[5]- Schelkunoff, S.
A., "Applied Mathematics for Engineers and Scientists", 2nd
Ed., p. 210, D. Van Nostrand Co., Princeton, New Jersey USA, 1965.
[6]- Roewe, D. and
Liou, K. N., J. Appl. Meteor., 17, 92, 1978.
[7]- Stamnes, K. S.,
Tsay, S. C., Wiscombe, W., and Jayaweera, K., App. Opt., 27, 2502, 1988.
