University of Southern California Dissertations and Theses

A matrix exponential approach to solving systems of linear differential equations with application to pharmacokinetic models

(USC Thesis Other)

A matrix exponential approach to solving systems of linear differential equations with application to pharmacokinetic models

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content A MATRIX EXPONENTIAL APPROACH TO SOLVING SYSTEMS
OF LINEAR DIFFERENTIAL EQUATIONS WITH APPLICATION
TO PHARMACOKINETIC MODELS
by
Gisela Spieler
A Thesis Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
MASTER OF SCIENCE
(Applied Mathematics)
December 1986
UMI Number: EP54448
All rights reserved
INFORMATION TO ALL USERS
The quality of this reproduction is dependent upon the quality of the copy submitted.
In the unlikely event that the author did not send a complete manuscript
and there are missing pages, these will be noted. Also, if material had to be removed,
a note will indicate the deletion.
UMI EP54448
Published by ProQuest LLC (2014). Copyright in the Dissertation held by the Author.
Microform Edition © ProQuest LLC.
All rights reserved. This work is protected against
unauthorized copying under Title 17, United States Code
ProQuest LLC.
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, Mi 48106- 1346
UNIVERSITY O F S O U T H E R N CA LIFO RN IA
THE GRADUATE SCH O O L
UNIVERSITY PARK
LO S A N G ELES. CA LIFO R N IA 9 0 0 0 7
This thesis, written by
Gisela Spieler
under the direction of h.^7:...Thesis Committee,
and approved by all its members, has been pre
sented to and accepted by the Dean of The
Graduate School, in partial fulfillment of the
requirements for the degree of
Master of Science
D ean
Date November 4, 1986
THESIS COMMITTEE
7 h a i? J h a n
WIDMUNG:
FÜR MEINEN VATER
"Danke, dap Du nie aufgegeben hast danach zu fragen."
11
A C K N O W L E D G E M E N T S
I would like to thank Professor Alan Schumitzky for his time and effort
as my advisor. My thanks to Professor David D'Argenio and Professor
Gary Rosen for their help and interest in the problem.
I would also like to thank the Department of Mathematics of the
University of Southern California for their financial support over the
past years.
Finally, a special thanks to my parents for all their encouragement and
support during my studies in the United States.
I l l
TABLE OF CONTENTS
Chapter
1. Introduction 1
2. Background 5
2.1 Definitions 5
2.2 Householder Transformations 5
2.3.1 QR - Factorization 6
2.3.2 QR - Algorithm 7
2.3.3 Shifts of Origin 9
2.4 Eigenvectors of Real Matrices by 12
QR Triangularizations
3. The Eispack Subroutines 15
3.1 ORTHES 15
3.2 HQR2 16
4. Application to Pharmacokinetic Systems 19
4.1.1 Linear Compartmental Models 19
4.1.2 Closed Systems
21
iv
4.1.3 Traps 21
4.1.4 Example 22
4.1.5 Eigenvalues of a Compartmental System 22
4.2 Application to Parameter Estimation 23
4.3 Diagram 25
5. Summary 26
References 27
Appendix 29
v
LIST OF FIGURES
Page
1. Figure 4.1 20
Diagram for a System with three Compartments
2. Figure 4.2 22
Example
3. Figure 4.3 25
Suggested Path for a Program solving IVP using ORTHES and
HQR2
VI
Chapter I
Introduction
The solution of a system of linear ordinary differential equations involves
the matrix exponential. In 1978 Cleve Moler and Charles Van Loan suggested
'T9 Dubious Ways" to calculate the matrix exponential (cf. [1]). This thesis
examines method #15 using an eigenvector approach. The eigenvectors are
computed via the QR algorithm.
More specifically, in many physical, biological and economic processes
we have systems of linear differential equations
x(t) = Ax(t)
where A is a given real or complex n x n matrix. A solution vector x(t) is
sought which satisfies the initial condition
x(0) = X q .
Theoretically the solution is given by
x(t) = e^bcQ,
where e^^ is the exponential of the matrix A, that can be defined as the
convergent power series
At T . . A 1 2 . 2
e — I + t A + A 4 - ... .
To calulate the matrix exponential, methods using matrix decompositions
proved to be efficient, especially for problems which require repeated
evaluation of e^^ These decompositions are based on similarity
transformations of the form
A = V D V k
The definition of the power series of e^‘ then implies that
e^t = Ve^W'^
The exponential of a diagonal matrix D = diag (1^,...,1^) is easily calculated
to be
e^t = diag (e^i\...,eV ).
Hence one possibility to calculate the matrix exponential is an eigenvalue
eigenvector approach. If 1^,...,1^are distinct eigenvalues of A, we can
calculate the eigenvectors v- corresponding to X ,-, i= l,.. .,n. Then Av- = X -v-
and if V = [Vj,.. .,v j we get
AV = VD,
where
D =diag (X^,...,XJ.
The approach considered in this thesis is especially attractive in conjunction
with parameter estimation problems.
Pharmacokinetic models are a mathematical description of kinetic phenomena
such as the absorption, distribution, metabolism and excretion of drugs.
Additionally these models are essential for designing individualized dosages to
achieve therapeutic goals. In most cases, to determine the model parameters
from input-output observations we have to solve a nonlinear parameter
estimation problem (cf. [14]).
The kinetics of a drug may be described by the following equation:
X = f(x,a,r,t)
x(0) = g(a)
y = h(x,a,r,t),
where x is the n-dimensional state vector, a the p dimensional time-invariant
parameter vector, r the k dimensional piece-wise constant input vector, y the
dimensional system output and t the independent variable. In general the
initial condition vector, x(0), can also be a function (g) of the unknown
parameter vector a. The parameter estimation problem is to find an estimate
of a which gives a "best fit" between the model and the measured data. For
example, one critérium for "best fit" could be the method of weighted least
squares. In this case the problem to be solved is, find the parameter vector a
which minimizes:
J(a) = % G . % w.. (y.(t, a) - z. (t.))^,
where y- and z- are the i-th elements of the model output and observation
vectors, respectively.
The term G-is the weight attached to the i-th output, while w-j are the
observation weights.
To solve for the least squares estimates we have to consider the following
problems:
(i) Integration of the model differential equations
(ii) Minimizing of the function J(a)
(iii) Analysis of the error associated with the parameter estimates and
subsequent model predictions.
An algorithm to solve the differential equations should be efficient since
numerous solutions of the model equations are generally required. Thus the
solution of linear time-invariant differential equations needs repeated
calculation of for many matrices A and time points t. In this thesis we
will outline a method to solve a linear system of ODE's using an eigenvector
approach. A brief description of the thesis now follows:
Chapter two wiU give the theoretical background for the transformation of
a general matrix to upper Hessenberg form and the QR algorithm. In chapter
three the Eispack subroutines HQR 2 and ORTHES are introduced. Chapter
four will present an application to pharmacokinetic modeling. Chapter five
summarizes advantages and problems of the algorithm. The appendix will
provide a listing of HQR2 and ORTHES.
Chapter 2
BACKGROUND
2.1 D efin ition s
i) A matrix H = (h-) e [R ™ is said to be upper Hessenberg if
hjj = 0 for i > j + 1.
ii) H is unreduced upper Hessenberg if all subdiagonal elements are non
zero.
iii) Q G is said to be orthogonal if Q^Q = I.
iv) Two matrices A and B are orthogonally similar if there exists an
orthogonal matrix Q such that Q^AQ = B.
v) An eigenvalue is defective if its algebraic multplicity exceeds its
geometric multiplicity.
2.2 H ouseholder transformations
A nxn matrix P of the form
P = I - 2vv^/yTy , for a nonzero vector v e IR ^ , is called a Householder
matrix.
Since P^ = I - 2w^/yTy and P^P = 1, P is symmetric and orthogonal.
In particular P can be constructed such that for any x e IR^, Px e span {e^}.
Since Px = x - 2 w^x/^T^ , choosing v = x + ae^ yields
v^x = x^x + axj and v^v = x^x+ 2ax^ + o?
and thus
T
2(x X + ax )
Px = x - — ----------------- (x + ae^) .
X X + 2aXj + a
For Px G sp{ej ) this implies a = ± I I x I I 2 and hence Px = ± I I x 1 1 2 6^.
Generalizing the above, we can use Householder matrices to "zero" any block
of vector components by choosing
2 2
a = x, +...+ X
k j
v'^ = [0,...,0, x^ + sgn(xj^) a,Xj^^p ... , Xj,0,...0].
Then
Px = [I - 2w"^/yTy] X = [Xp...,Xj^_p - sgn(Xj^)a,0,..0,Xj^j, ..., x j
and P can be written as
P = diag P,
where P is a Householder matrix generated by v's nontrivial part.
2.3.1 OR - factorization
The QR algorithm is based on the fact that each nxn matrix B can be
decomposed into B = QR, with unitary matrix Q and upper triangular matrix
R.
We can calculate Q by applying a finite number of Householder matrices to B
such that
^n-i • •••Pj • Pq • B = Q^B = R.
Assume after r steps that B^ has the following blockstmcture:
B =
r
u
0
11
u
12
u
22
with Ujj an rxr upper triangular matrix.
Defining
Pr =
r,r
0
we obtain
11
0
I -2v V V v
n-r,n-r r r r r
0
Hence V mustbe chosen such Aat the first column o f f 1 -2v v V v )U„„ has all zero entries except
r \ n -r ,n -r x x x x/ L L ^
the first. After (n-1) steps, B ^ ^ j = P^ ........Pq B is triangular and therefore
QB^^ = B is the required factorization.
2.3.2. OR Algorithm
The QR-Algorithm reduces a matrix to quasi-tiiangular form by orthogonal
similarity transformation.
The theorem on the real Schur decomposition of a matrix A g IR ^ (cf. [2])
provides us with the existence of an orthogonal matrix Q g IR ^ such that
Q^AQ =
U u .
11« Im
0 u
where each is a 1x1 matrix or a 2x2 matrix with conjugate complex
eigenvalues.
To obtain Q from A we form an infinite sequence of orthogonal matrices Q j^
and upper triangular matrices R j^ . defined by the following relations :
Aq = A for k=0,l,...
= Qk %
^k+l = % Qk
■ n'usA ^^i=R ,Q k = Q ^Q ,R ,Q fc
= Q k \ Q k -
Each Ay. is similar to Aq = A. Qq can be chosen arbitrarily.
For each iteration a QR - factorization of A j^ has to be performed. The amount
of calculations causes the QR algorithm to be impractical, unless the matrix A
has many zero entries.
For example, if A is of upper Hessenberg form the algorithm becomes useful.
Since Qq is arbitrary, it is reasonable to choose our orthogonal Qq such that
Q J a Qq = H is of upper Hessenberg form.
Qq is obtained from a finite product of Householder matrices in a similar way
as in the QR factorization, (cf. [2]).
Multiplication of an upper Hessenberg matrix with a triangular matrix
preserves the Hessenberg structure, and hence all = R j^ Q j^
= % Ak RiT
are of upper Hessenberg form.
The proof of the convergence of {A ^^} to upper triangular form can be found
in [2],[4].
2.3.3 Shifts of Origin
The following relations describe Francis's QR Algorithm with shifts of origin
(c /. [9])
Qk %
^k+1 ^ ^k Qk +
Thus A jj. and are orthogonally similar:
\ + i =BkQk + V
= Q l \ Q k
Performing two steps of the QR algorithm at the same time supplies an
economical way to keep the arithmetic real.
After two steps we have
^k+2 " Qk+1 Qk "^k Qk ^k+l
giving
Qk Qk+l ^k+2 = ^ k Qk Qk+1
Qk Qk+, Rk+i \ = Qk ( I ) Rk
=Qk(RkQk + V - W ) R k
“ Qk Rk (Qk Rk V ■ \ + l O
=(Ak - V )(A k -
Defining Q = QkQk+i>
R = Rk+iRk“ d
n = (Ak - V )(A k - V i D .
we obtain
n=QR
and
Ak+2 = Q^AkQ.
Q is a real matrix since 7^ and are chosen both real or complex conjugate
and hence the product H is real. The shifts of origin at each step are taken to
10
be the two roots of the 2x2 matrix in the right hand comer of the current
where
a(k)
n -l,n -2 n-1,n-1 n -l,n
n,n
The shifts satisfy
=
+1 n -l,n -l
+ a
(k)
n,n
(k)
. a
(k)
\ ' ^ V l,n - 1 n,n
With these shifts the convergence of the subdiagonal elements is rapid and at a
certain stage the
(k) (k) ....
\ n - i o r are negligible.
(k) (k)
If dn n 1 is very small, then a^^ can be considered an eigenv alue of the original matrix and the matrix
deflates. That is the next step can be performed with the smaller Hessenberg
matrix obtained by dropping the last row and column.
11
(k)
If 1 n 2 is sm all, th e tw o eigenvalues of the 2x2 principal subm atrix in theright hand com er are taken
to be eigenvalues of A and we can drop the last two rows and columns.
After each two step iteration all subdiagonal elements are examined to see if
any of them are negligible and the Hessenberg matrix decouples into smaller
submatrices.
A sufficient condition for A ^ ^ to converge to quasi-triangular form is that A has
n distinct eigenvalues such that \Xy I > I X 2 I >...> I I > 0 (cf.[4]).
It can be shown (cf. [4],[5]) that the a^ • > 0 approximately like (X/X-
Thus with a shift of
k
(
X - g \
— — I tends to zero rapidly.
V rti J
2.4 Eigenvectors of Real Matrices bv OR Triangularizations
Given an upper Hessenberg matrix H and a quasi-triangular matrix T such that
Q^HQ = T. Then the eigenvalues of H are the eigenvalues of the 1x1 and 2x2
principal submatrices of T.
If a block on the diagonal of T has a pair of real eigenvalues in row and
column r and r+1, the true triangular form can be established by a similarity
transformation involving a real rotation in the (r,r 4-1) plane.
12
For a real eigenvalue X -, the corresponding eigenvector x- may be calculated as
follows
(1) x.j = 0 j = i+ l,...,n
X ii = 1
-ij = - ( i
j = i-l,...l.
If X - = X j (j < i) the divisor is replaced by an expression which ensures that
we obtain the requisite number of independent eigenvectors. For this case
and for eigenvalues corresponding to non-linear divisors refer to [7].
A 2x2 block with complex conjugate eigenvalues can not be triangularized
without complex arithmetic.(Note: most compartmental models have only real
eigenvalues, see Chapter 4.)
However, the eigenvectors can still be calculated by backsubstitution. It is
sufficient to determine an eigenvector corresponding to the eigenvalue with
positive imaginary part, since the other vector is just the complex conjugate.
The components X j^ and X j ^. corresponding to a 2x2 block in rows j and j-1
are given by
13
1
" " " ^ ^ j-1,A i
k=j+l
i
S o -l^ j-l,i ^ ^jk^ki
k=j+l
Taking x^ = 1 -I- Oi, the component x^^ - satisfies
(^i-i,i-i “ ^i) ^i-i,i + Vi,i ^
^i,i-l ^i-l,i (^ii ■ ^i) -
The remaining components are given by (1) and (2).
It is worth noting that the matrix R containing the eigenvectors obtained above
is upper triangular.
14
Chapter 3
The Eispack Subroutines
3.1 ORTHES
Reduction of a real general matrix to upper Hessenberg form using orthogonal
similarity tranformations.
Given a matrix A of the form
A =
where T,R are upper triangular matrices and B is a pxp square matrix in rows
e through i. The subroutine ORTHES reduces B to upper Hessenberg form.
Assume B is already of upper Hessenberg form in columns 1 through r-1. To
avoid possible underflow the elements below the diagonal in column r are
scaled.
T X Y
O B Z
O O R
Define
,2 . .2
= ^+l,r + ---+ d p .
u =
1/2
= [ 0 0 'b r+ l.r (b r+l.r) V
H = a - sgn a
1/2
15
Then, the Householder matrix = I - n^n^/H^ is orthogonal and symmetric.
The multiplication
Pj. Pj. does not effect the zeros in columns one through r-1. B ^.^^ = P^ B^
P^ is computed via
C r+l = P r B , =
B r+1 = Qf+l P r = Q r + 1 ' (^,+1 "P/H, U^.
The above steps are repeated until B is of upper Hessenberg form. The
information about the orthogonal transformations is stored in the vector ORT
and the lower triangle of A. The subroutine ORTRAN accumulates these
orthogonal transformations into the array Z.
3 .2 HOR2
The subroutine HQR2 calculates eigenvalues and eigenvectors of a real
upper Hessenberg matrix H using the QR method. A sequence of upper
Hessenberg matrices is formed iteratively. Each element of the sequence
which converges to quasi-triagular form T is orthogonally similar to the
original matrix H.
The eigenvalues of H are then equal to the eigenvalues of the 1x1 and 2x2
principal submatrices of T. Before each iteration the matrix is checked for
possible deflation that would be partitioning of the matrix into smaller blocks.
The algorithm can be discribed as follows:
16
Repeat:
Depending on the relative machine precision e, all subdiagonal elements
satisfying
I a - I < 6 ( I a. J + I a - _ j I ) are set equal to zero.
The smallest non negative number p and the largest non negative number q are
determined such that H can be described as
H 12 Hi3
P
H =
0
^ 2 2 ^23
n-p-q ,
0 0
^33
q
P
n-p-q
q
where H22 is an unreduced upper Hessenberg matrix and H33 is upper
quasi-triangular.
If q = n (ie H is quasi-triangular) triangularize all 2x2 diagonal blocks with
real eigenvalues, accumulate the transformations and quit.
Otherwise apply the Franics QR step with double shift of origin to the
unreduced Hessenberg matrix H22. The shifts of origin a^ and a2 are chosen
to be the real or complex conjugate eigenvalues of the tailing 2 x2 principal
submatrix of H22.
Determine a product of Householder matrices Z such that Z^(H22 - a^I)(H22 -
a2l) is upper triangular and overwrite H22 by
17
Once the lowest 2x2 matrix splits, its eigenvalues are taken to be the
eigenvalues of H and we go back to repeat above until all eigenvalues are
found.
HQR2 then calculates the eigenvectors of the quasi-triangular matrix T by
back substitution and transforms them to eigenvectors of the original matrix H
using the information about the orthogonal similarity transformation.
HQRORT, the modified version of HQR2, completes the orthogonal
reduction of H to quasi-triangular form and accumulates the remaining
transformations. The statement labeled 60 in HQR2 is changed to 60 If (en
.It. LOW) Go to 1001 (cf. [3]).
18
Chapter 4
Application to Pharmacokinetic Systems
4.1.1 Linear Compartmental Models
An n-th order linear system
(1) x(t) = Ax(t) + b(t)
x(0 ) = X q t > 0
with constant coefficient matrix A e , such that
a) each diagonal element is nonpositive
b) every off diagonal element is nonnegative
c) the sum of column j is the non-positive number - a^j, j= l,.. .,n
is called a n-compartmental system with state (or compartmental) matrix A,
state variable x(t) = [x^(t),... ,x^(t)]^ and
input vector (forcing function) b(t) = [b^(t),...,b^(t)]^
Figure 4.1
Diagram for a system with three compartments, where k-j is the rate of
exchange from compartment j into compartment i, i^}.
19
b,(t)
b,(t)
The state matrix A for such a linear time invariant compartmental system is
A =
In general: a-- = k-- i j; i,j = l,...,n and
a..
i=i
20
4.1.2 C losed System s
If no compartment excretes material to the environment, (1) is said to be a
closed system. In a closed system Icq j = 0, V j. This implies that all column
sums are equal to zero and zero is an eigenvalue of A, (cf. [10]).
4.1.3 T ra p s
The system contains a trap (sink) if there is a subset S of compartments such
that
Vj e S, V i ^ S =» a- = 0,
that is, a trap is a subsystem with no excretion to things outside itself,
(cf. [10]).
Proposition: A system containing a trap has a zero eigenvalue.
P roof:
Assume WLOG that compartments 1 through k form a trap of the system, i.e.
for all je {l,...,k} and for all i e {k+1,. ..,n},
a-j = 0 .
Hence the state matrix can be partitioned
A
"^11 ^ 1 2
0 A
22
n-k
n-k
21
and det A = det A^^# det A22. But, A^^ is the compartment matrix of a
closed system and thus has a zero eigenvalue and det A^j = 0 => det A = 0.
4.1.4 E xam ple
Consider a closed system with 4 compartments and a trap.
"21
V
k i3 ^ 1 ^24 ^ 42
>
>
k
43
V
3 4
The compartments 2 and 4 are forming a trap and X = 0 is eigenvalue of
multiplicity one.
4.1.5 Eigenvalues of a compartmental system
A sufficient condition for the eigenvalues of a compartmental system to be real
is that all cycles of r compartments within a system (a cycle being a series of
interconnected compartments
ap'
22
Some of the eigenvalues can be complex but, in the linear time invariant
model, any oscillation is heavily damped. For this reason, the eigenvalues are
usually negative real except for closed models or models containing a sink
which have an eigenvalue at the origin (cf. [13]).
4.2 Application to Parameter Estimation
Parameter estimation in pharmacokinetic modeling gives rise to systems of
linear ordinary differential equations
x(t) = x(t)
(IVP)
X(0) = X q
t > 0
for a family ( Aj.} of compartmental matrices. This system can be solved in
the following way.
i) Reduce A to upper Hessenberg form by orthogonal similarity
transformations
Q/'AQj =H .
ii) Find an orthogonal matrix Q2 and a quasi-triangular matrix T such that
Q2 HQ2 = T
iii) Calculate the eigenvectors of T by backsubstitution to obtain a matrix
R satisfying
TR = RD, where D is diagonal matrix with the eigenvalues of A on its
diagonal. R can be taken to be a triangular matrix. 23
Hence the matrix Z containing the eigenvectors of A is Q«R, where Q = Q1Q2
is orthogonal, iv) Solve Ry^ = Q^Xq for yQ by backsubstitution.
Thus the initial value problem has the solution x(t) = QRe^^y^.
Once Q, R and y^ are computed, the solution x(t) can be evaluated for
different times t- with relatively little work.
To estimate the accuracy of y^, we have to consider cond(QR) = cond(R) in
the 2-norm. If cond(R) is too large this algorithm wiU not solve the problem.
24
4.3 Suggested path for a program solving IVP using ORTHES and HQR2
Read A, read
read eps
c jl ORTHB & CRTRAN
Q^AQ = H
caüHQR2
TR = RD
ccnd(R) < ecB
Solve Ryg = ^Xg
quit
D u
x(t.)=ZRe 7 (
1 = i + 1
yes
Ifi< k
-2-5J
Chapter 5
Summary
The presented algorithm to calculate the matrix exponential has the advantage
that most of the calculations are done independently of t, while for example
W ard’ s method using Fade approximations (c /. [8]), has to perform a new
exponential calculation for each value of t even though the matrix A is
unchanged.
The eigenvectors calculated in HQR2 corresponding to nondefective multiple
eigenvalues are linearly independent. If a multiple eigenvalue is defective,
HQR2 produces linearly dependent eigenvectors. This will reflect in the
cond(R) and indicate that the algorithm can not solve the problem. However
in our pharmacokinetic applications, the eigenvalues are in general distinct,
except for multiple zeros which were in most examples nondefective.
26
References
[1] C Moler, C. van Loan: 19 Dubious Ways to Calculate the Exponential of
a Matrix, SIAM Review vol. 20 No. 4 (1978)
[2] G.H. Golub, C van Loan: Matrix Computations, Johns Hopkins Univ.
Press, Baltimore 1983
[3] B.T. Smith: Matrix Eigensystem Routines; EISPACK guide, 2nd ed.
Lecture notes in Comp Science, Springer-Verlag, NY 1976
[4] E.K. Blum: Numerical Analysis and Computation Theory and Practice,
Addison Wesley 1972
[5] B. Parlett: Global Convergence of the Basic QR Algorithm on
Hessenberg Matrices, Math of Comp 22, pp. 803-817 (1968)
[6] S. M artin, G. Peters, J.H. W ilkinson: QR Algorithm for Real
Hessenberg Matrices, Num Math 14, pp. 219-231 (1970)
[7] G. Peters, J.H. Wilkinson: Eigenvectors of Real and Complex Matrices
by LR and QR Triangularizations, Num Math 16, pp. 181-204 (1970)
[8] R.C. Ward: Numerical Computation of the Matrix Exponential with
Accuracy Estimate, SIAM J. Num Anal, 14, pp. 600-610 (1977).
[9] J.G.F. Francis: The QR Tranformation, I, II, Comp J. V.4 1960/1961
pp. 265-271, 332-345.
[10] D.H. Anderson: Compartmental Modeling and Tracer Kinetics, Springer
Verlag Berlin 1983.
[11]U. Godfrey: Compartmental Models and Their Application, Academic
Press, NY, 1983. 27
[12] J.A. Jacques: Compartmental Analysis in Biology and Medicine,
Elsevier Amsterdam 1972.
[13] R.F. Brown: Compartmental System Analysis: State of the Art, IEEE
Trans on Biomed Eng. Vol. BME 27 No. 1 pp. 1 - 11 (1980).
[14] D. DArgenio, A. Schumitzky: A Program Package for Simulation and
Parameter Estimation in Pharmacokinetic Systems, Comp. Programs in
Biomed. 9(1979) pp. 115-134.
28
A ppendix
subroutine orthes(nm,n,low,igh,a,ort)
c
integer i,j,m,n,ii,Jj,la,mp,nm,igh,kpl,low
double precision a(nm,n),ort(igh)
double precision f,g,h,scale
c
c this subroutine is a translation of the algol procedure orthes,
c n'um. math. 12, 349-358(1958) by martin and wilkinson.
c handbook for auto, comp., vol.ii-linear algebra, 339-358(1971).
c
c given a real general matrix, this subroutine
c reduces a submatrix situated in rows and columns
c low through igh to upper hessenberg form by
c orthogonal similarity transformations,
c
c on input
c
c nm must be set to the row dimension of two-dimensional
c array parameters as declared in the calling program
c dimension statement.
c
c n is the order of the matrix,
c
c low and igh are integers determined by the balancing
c subroutine balanc. if balanc has not been used,
c set low=l, igh=n.
c
c a contains the input matrix,
c
c on output
c
c a contains the hessenberg matrix. information about
c the orthogonal transformations used in the reduction
c is stored in the remaining triangle under the
c hessenberg matrix,
c
c ort contains further information about the transformations,
c only elements low through igh are used,
c
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
this version dated august 1983
la = igh - 1
kpl = low + 1
if (la .It. kpl) go to 200
do 180 m = 1^1, la
h = 0.OdO
ort (m) = 0 . OdO
scale = 0.OdO
......... scale column (algol toi then not needed)
do 90 i = m, igh
90 scale = scale + dabs (a (i,m-l))
if (scale .eg. 0.OdO) go to 180
mp = m + igh
......... for i=igh step -1 until m do -- .....
do 10 0 ii = m, igh
i = mp - ii
ort (i) = a (i,m-l) / scale
h = h + ort (i) * ort (i)
100 continue
29
c
g = -dsign (dsqrt (h) , ort (m) )
h = h - ort (m) * g
ort (m) = ort (m) - g
.......... form (i-(u*ut)/h) * a .........
do 130 i = m, n
f = O.OdO
......... for i=igh step -1 until m do --
do 110 ii = m, igh
i = mp - ii
f = f + ort (i) * a (i, ])
110 continue
f = f / h
do 120 i = m, igh
120 a (i, i) = a(i,j) - f * ort (i)
130 continue
......... form (i- (u*ut) /h) *a* (i- (u*ut) /h)
do 160 i = 1, igh
f = O.OdO
for j=igh step -1 until m do --
do 140 ]] = m, igh
j = mp - ii
f = f + ort (i) * a(i,i)
140 continue
c
f = f / h
do 150 i = m, igh
150 a (i, i) = a (i, i), - f * ort (i)
c
160 continue
c
ort (m) = scale * ort (m)
a(m,m-l) = scale * g
180 continue
c
200 return
end
30
subroutine hqrZ (nm, n, low, igh, h, wr, wi, z, ierr)
c
integer i,j,k,l,m,n,en,ii,jj,ll,mm,na,nm,nn,
X igh,itn,its,low,mp2,enm2,ierr
double precision h (nm, n) , wr (n) , wi (n) , z (nm, n)
double precision p, q, r, s, t, w, x,y, ra, sa, vi, vr, zz,norm, tstl, tst2
logical notlas
c
c this subroutine is a translation of the algol procedure hqr2,
c num. math. 16, 181-204(1970) by peters and wilkinson.
c handbook for auto, comp., vol.ii-linear algebra, 372-395(1971).
c
c this subroutine finds the eigenvalues and eigenvectors
c of a real upper hessenberg matrix by the qr method. the
c eigenvectors of a real general matrix can also be found
c if elmhes and eltran or orthes and ortran have
c been used to reduce this general matrix to hessenberg form
c and to accumulate the similarity transformations.
c
c on input
c
c nm must be set to the row dimension of two-dimensional
c array parameters as declared in the calling program
c dimension statement,
c
c n is the order of the matrix,
c
c low and igh are integers determined by the balancing
c subroutine balanc. if balanc has not been used,
c set low=l, igh=n.
c
c h contains the upper hessenberg matrix,
c
c z contains the transformation matrix produced by eltran
c after the reduction by elmhes, or by ortran after the
c reduction by orthes, if performed. if the eigenvectors
c of the hessenberg matrix are desired, z must contain the
c identity matrix,
c
c on output
c
c h has been destroyed,
c
c wr and wi contain the real and imaginary parts,
c respectively, of the eigenvalues. the eigenvalues
c are unordered except that complex conjugate pairs
c of values appear consecutively with the eigenvalue
c having the positive imaginary part first. if an
c error exit is made, the eigenvalues should be correct
c for indices ierr+1,...,n .
c
c z contains the real and imaginary parts of the eigenvectors,
c if the i-th eigenvalue is real, the i-th column of z
c contains its eigenvector. if the i-th eigenvalue is complex
c with positive imaginary part, the i-th and (i+1)-th
c columns of z contain the real and imaginary parts of its
c eigenvector. the eigenvectors are unnormalized. if an
c error exit is made, none of the eigenvectors has been found,
c
c ierr is set to
c zero for normal return,
c j if the limit of 30*n iterations is exhausted
c while the ]-th eigenvalue is being sought,
c
c calls cdiv for complex division,
c
31
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
c this version dated august 1983.
c
ierr = 0
norm = 0.OdO
k = 1
c ......... store roots isolated by balanc
c and compute matrix norm .
do 50 i = 1, n
c
do 40 j = k, n
40 norm = norm + dabs(h(i,j))
c
k = i
if (i .ge. low .and. i .le. igh) go to 50
wr (i) = h(i, i)
wi (i) = O.OdO
50 continue
c
en = igh
t = 0.OdO
itn - 30*n
c ......... search for next eigenvalues .........
50 if (en .It. low) go to 340
its = 0
na = en - 1
enm2 = na - 1
c ......... look for single small sub-diagonal element
c for l=en step -1 until low do -- .....
70 do 80 11 = low, en
1 = en + low - 11
if (1 .eq. low) go to 10 0
s = dabs (h(l-l, 1-1) ) + dabs (h(l, 1) )
if (s .eq. O.OdO) s = norm
tstl = s
tst2 = tstl + dabs (h(l, 1-1) )
if (tst2 .eq. tstl) go to 100
80 continue
c ......... form shift ..........
100 X = h (en,en)
if (1 .eq. en) go to 270
y = h(na,na)
w = h (en, na) * h(na.en)
if (1 .eq. na) go to 280
if (itn .eq. 0) go to 1000
if (its .ne. 10 .and. its .ne. 20) go to 130
c ......... form exceptional shift .........
t = t + X
c
do 120 i = low, en
120 h(i,i) = h(i,i) - x
c
s = dabs (h (en,na) ) + dabs (h (na, enm2) )
X = 0.75dO * s
y = X
w = -0.4375dO * s * s
130 its = its + 1
itn = itn - 1
c ......... look for two consecutive small
c sub-diagonal elements.
c for m=en-2 step -1 until 1 do -- ..........
do 140 mm = 1, enm2
32
m = enm2 + 1 - mm
zz = li(m,m)
r = X - zz
s = y - zz
p = (r * s - w) / h(m+l,m) + h(m,m+l)
q = h(m+l,m+l) - zz - r - s
r = h(m+2,m+l)
s = dabs (p) + dabs (q) + dabs (r)
P = P / s
q = q / s
r = r / s
if (m .eq. 1) go to 150
tstl = dabs (p) * (dabs (b(m-l,m-l) ) + dabs (zz) + dabs (h (m+l,m+l) ) )
tst2 = tstl + dabs (h (m,m-l) ) * (dabs (q) + dabs (r) )
if (tst2 .eq. tstl) go to 150
140 continue
c
150 mp2 = m + 2
c
do 150 i = np2, en
h(i, i-2) = 0 .OdO
if (i .eq. mp2) go to 160
h(i,i-3) = O.OdO
160 continue
c..... .......... double qr step involving rows 1 to en and
c columns m to en ..........
do 260 k = m, na
notlas = k .ne. na
if (k .eq. m) go to 170
p = h(k,k-l)
q = b(k+l,k-l)
r = O.OdO
if (notlas) r = h(k+2,k-l)
x = dabs (p) + dabs (q) + dabs (r)
if (x .eq. O.OdO) go to 260
p = p / X
q = q / X
r = r / X
170 s = dsign(dsqrt (p*p+q*q+r*r) ,p)
if (k .eq. m) go to 180
h(k,k-l) = -s * X
go to 190
180 if (1 .ne. m) h(k,k-l) = -h(k,k-l)
190 p = p + s
X = p / s
y = q / s
zz = r / s
q = q / P
r = r / p
if (notlas) go to 225
c..... ......... row modification ....
do 200 j = k, n
p = b(k, j) + q * h(k+l, j)
b(k,j) = b(k, ]) - p * X
h (k+1, ]) = b(k+l ,j) - p * y
200 continue
c
j = min0(en,k+3)
c.............. column modification . .......
do 210 i = 1, j
p = X * h(i,k) + y * h(i,k+l)
h(i,k) = h(i,k) - p
b(i,k+l) = h(i,k+l) - p * q
210 continue
c..... ......... accumulate transformations .........
do 220 i - low, igh
33
p=x*z(i,k) + y* z(i,k+l)
z (i,k) = z (i,k) - p
z(i,k+l) = z(i,k+l) - p * q
220 continue
go to 255
225 continue
c..... .......... row modification ..........
do 230 i = k, n
p = h(k,j) + q * h(k+l,j) + r * h(k+2,J)
h(k, j) = h(k, j) - p * X
h(k+l, ]) = h(k+l,j) - p * y
h(k+2,j) = h(k+2, j) - p * zz
230 continue
c
i = minO(en,k+3)
c ......... column modification .........
do 240 i = 1, ]
p = X * h (i,k) + y * h (i, k+1) + zz * h(i,k+2)
h(i,k) = h(i,k) - p
h(i,k+l) = h(i,k+l) - p * q
h ri,k+2) = h(i,k+2) - p * r
240 continue
c .......... accumulate transformations .......
do 250 i = low, igh
p = X * z(i,k) + y * z (i,k+1) + zz * z(i,k+2)
z (i,k) = z (i,k) - p
z (i,k+l) = z (i,k+l) - p * q
z(±,'k+2) = z(i,k+2) - p * r
250 continue
255 continue
c
26G continue
c
go to 70
c ......... one root found .........
270 h(en,en) = x + t
wr (en) = h (en, en)
wi(en) = 0.OdO
en = na
go to 60
c ......... two roots found .........
280 p = (y - x) / 2.OdO
q = p * p + w
zz = dsqrt (dabs (q) )
h (eh, en) = x + t
X - h (en, en)
h (na, na) = y + t
if (q .It. O.OdO) go to 320
c.............. real pair .........
zz = p + dsign(zz,p)
wr (na) = x + zz
wr (en) = wr (na)
if (zz .ne. O.OdO) wr (en) = x - w / zz
wi (na) = 0 . OdO
wi (en) = 0 . OdO
X = h (en, na)
s = dabs (x) + dabs (zz)
p = X / s
q = zz / s
r - dsqrt (p*p+q*q)
p = p / r
q = q / r
c ......... row modification .........
do 290 i - na, n
zz = h(na, ])
h(na, i) = q * zz + p * h(en,j)
3ÀJ
h(en,j) = q * h (en, j) - p * zz
290 continue
c
do 300 i = 1, en
zz = h(i,na)
h (i, na) = q * zz + p * h(i, en)
h(i,en) = q * h (i, en) - p * zz
300 continue
c ......... accumulate transformations .........
do 310 i = low, igh
zz = z(i,na)
z (i, na) = q * zz + p * z (i, en)
z (l,en) = q * z (i, en) - p * zz
310 continue
c
go to 330
c ....... complex pair .........
320 wr (na) = x + p
wr (en) = x + p
wi (na) = zz
wi (en) - -zz
330 en = enm2
go to 60
c ......... all roots found, backsubstitute to find
c vectors of upper triangular form .........
340 if (norm .eq. O.OdO) go to 1001
c ......... for en=n step -1 until 1 do -- .........
do 800 nn = 1, n
en = n + 1 - nn
p = wr (en)
q = wi (en)
na = en - 1
if (q) 710, 600., 800
c . ......... real vector .........
600 m = en
h (en,en) = 1.OdO
if (na .eq. 0) go to 800
c ......... for i=en-l step -1 until 1 do -- ..........
do 700 ii = 1, na
i = en - ii
w = h(i, i) - p
r = 0.OdO
c
do 510 j = m, en
610 r = r + h(i,j) * h(i, en)
c
if (wi (i) .ge. O.OdO) go to 630
zz = w
s = r
go to 700
630 m = i
if (wi (i) .ne. O.OdO) go to 640
t = w
if (t .ne. O.OdO) go to 635
tstl = norm
t = tstl
632 t = 0.OldO * t
tst2 = norm + t
if (tst2 .gt. tstl) go to 632
635 h (i, en) = -r / t
go to 680
c
640 X = h(i, i+1)
y = hri+1, i)
q = (wr (i) - p) * (wr (i) - p) + wi (i) * wi (i)
t = (x * s - zz * r) / q
3 5
h(i,en) = t
if (dabs (x) .le. dabs (zz) ) go to 650
h (i+1, en) = (-r - w * t) / x
go to 680
650 h (i+1,en) = (-s - y * t) / zz
c
c ......... overflow control .........
680 t = dabs (b(i, en) )
if (t .eq. O.OdO) go to 70 0
tstl = t
tst2 = tstl + 1.OdO/tstl
if (tst2 .gt. tstl) go to 700
do 690 j = i, en
h(j,en) = h(j,en)/t
690 continue
c
700 continue
c end real vector ;.....
go to 800
c ......... complex vector .......
710 m = na
c ......... last vector component chosen imaginary so that
c eigenvector matrix is triangular .........
if (dabs (h (en, na) ) .le. dabs (h (na, en) ) ) go to 720
h(na,na) = q / h (en, na)
h (na, en) = - (h(en, en) - p) / h (en,na)
go to 730
720 call cdiv (0 . OdO ,-h (na, en) ,h (na, na)-p, q,h (na,na) ,h(na, en) )
730 h (en, na) = O.OdO
h (en, en) = 1. OdO
enm2 = na - 1
if (enm2 .eq. 0) go to 800
c ......... for i=en-2 step -1 until 1 do -- . •........
do 795 ii = 1, enm2
i = na - ii
w = h(i,i) - p
ra = 0.OdO
sa = 0.OdO
c
do 760 j = m, en
ra = ra + h(i,j) * h(j,na)
sa = sa + h (1, i) * h(j,en)
760 continue
if (wi (i) .ge. O.OdO) go to 770
zz = w
r = ra
s = sa
go to 795
770 m = i
if (wi (i) .ne. O.OdO) go to 780
call cdiv (-ra, -sa, w, q,h (i, na) ,h(i, en) )
go to 790
.......... solve complex equations . . . .....
780 X = h(i,i+l)
y = h(i+l, i)
vr = (wr (i) - p) * (wr (i) - p) + wi (i) * wi (i) - q
vi = (wr (i) - p) * 2. OdO * q
if (vr .ne. O.OdO .or. vi .ne. O.OdO) go to 784
tstl - norm * (dabs (w) + dabs (q) + dabs (x)
X + dabs (y) + dabs (zz) )
vr = tstl
783 vr = 0.OldO * vr
tst2 = tstl + vr
if (tst2 .gt. tstl) go to 783
784 call cdiv(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra, vr, vi.
36
X h (i, na) ,h (i, en) )
if (dabs (x) .le. dabs (zz) + dabs (q) ) go to 785
h (i+l,na) = (-ra - w * h(i,na) + q * h (i, en) ) / x
h (i+1, en) = (-sa - w * hri,en) - q * h(i,na) ) / x
go to 790
785 call cdiv (-r-y*h(i,na) ,-s-y*h(i,en) ,zz,q,
X h (i+l,na) ,h (i+1, en) )
c
c ......... overflow control .........
790 t = dmaxl (dabs (h (i,na) ) , dabs (h (i, en) ) )
if (t .eq. O.OdO) go to 795
tstl = t
tst2 = tstl + 1.OdO/tstl
if (tst2 .gt. tstl) go to 795
do 792 j = i, en
h(j,na) = h(j,na)/t
b(j,en) = h(j,en)/t
792 continue
c
795 continue
c end complex vector ........
800 continue
c ......... end back substitution.
c vectors of isolated roots .......
do 840 i = 1, n
if (i .ge. low .and. i .le. igh) go to 840
c
do 820 j = i, n
820 z (i, j) = h(i, j)
c
840 continue
c ......... multiply by transformation matrix to give
c vectors of original full matrix.
c for j=n step -1 Until low do -- ..........
do 880 jj = low, n
j = n + low - j j
m = minO (] , igh)
c
do 880 i = low, igh
zz = O.OdO
c
do 850 k = low, m
860 zz = zz + z(i,k) * h(k,j)
c
Z(i,j) = zz
880 continue
c
go to 1001
c ......... set error -- all eigenvalues have not
c converged after 30*n iterations .........
1000 ierr = en
1001 return
end
- 17_

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Semi-linear regression problems with an application in pharmacokinetics

PDF

On the multiple model adaptive control of pharmacokinetic systems

PDF

The extended Kalman filter as a parameter estimator with application to a pharmacokinetic example

PDF

The application of an approximately optimal nonlinear stochastic control and estimation algorithm to a pharmacokinetic problem

PDF

The NONMEM approach to estimation of population pharmacokinetic parameters

PDF

Estimation of states and parameters in nonlinear pharmacokinetic models by extended Kalman filtering

PDF

Estimating the lifetime risk of a late onset disease when sampling sibships with probability proportional to size: Testing autosomal dominant and recessive genetic hypotheses

PDF

An application of the ridge regression method to a neurological nursing care study

PDF

The asymptotic behavior of solutions of a second order linear differential equation

PDF

An algorithm for solving differential delay equations

PDF

A matrix algebra approach for determining the stationary random responses of linear structural systems

PDF

A comparison of findings from City of Hope's Mobile Screening Unit program with data from the Los Angeles County Cancer Surveillance Program

PDF

Stabilization of a linear system via rotational control

PDF

A course of study in orthopedic nursing

Systematic biochemical analysis of single nucleotide polymorphisms in the 3<super>'</super> untranslated region of the SRD5A2 gene in the context of prostate cancer

PDF

Systematic biochemical analysis of single nucleotide polymorphisms in the 3' untranslated region of the SRD5A2 gene in the context of prostate cancer

PDF

An evaluation of a concept formation test for two and three year old deaf and hard of hearing pupils

PDF

Nursing management of elderly patients with delirium and/or dementia

PDF

Herpesvirus and HIV type 1 infections: A seroepidemiologic study of relationships in selected populations

PDF

Innervation of the human tooth

PDF

A case-control study of maternal risk factors for thyroid cancer in young women

Asset Metadata

Creator Spieler, Gisela (author)

Core Title A matrix exponential approach to solving systems of linear differential equations with application to pharmacokinetic models

Contributor Digitized by ProQuest (provenance)

Degree Master of Science

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag applied sciences,Health and Environmental Sciences,OAI-PMH Harvest

Format application/pdf (imt)

Language English

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c37-173140

Unique identifier UC11647960

Identifier EP54448.pdf (filename),usctheses-c37-173140 (legacy record id)

Legacy Identifier EP54448.pdf

Dmrecord 173140

Document Type Thesis

Format application/pdf (imt)

Rights Spieler, Gisela

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA