Annals of Biomedical Engineering, Vol. 18, pp. 597-621, 1990 Printed in the USA.All rights reserved.
0090-6964/90 $3.00 + .00 Copyright9 1990PergamonPress plc
Fractal System-A Time Domain Approach H. H. Sun and A. Charef Department of Electrical and Computer Engineering Drexel University Philadelphia, PA (Received 3/9/90; Revised 5/24/90)
A method to analyze thefractal system in the time domain is presented so that the dynamic behavior o f the system can be studied. The fractal system is represented by a set o f linear time-varying differential equations whose order depends on the order o f the system under non-fractal condition. Four different types o f fractal system are considered and their solutions in the time domain are presented. These analyses show that the fraetal system is dynamically more stable with smooth changes o f magnitude and less oscillatory than the non-fractal system. Examples o f the physiological system o f the conduction path ways in the heart and also the polarization phenomena o f noble metal are presented to illustrate the phonomena. Keywords-Fractal systems, Bioelectrode, Cardiovascular system.
INTRODUCTION Numerous studies have been made about the fractal system ever since the phenomenon of 1/f-noise spectrum was first introduced by Van Der Ziel back in 1950 (20). During the past decade we have seen a tremendous growth of study in this subject especially toward the applications in meteorology, astronomy, magnetization, electrochemistry, geology, physiology, and many other natural phenomena (12). Most of these studies consider the fractal system as a r a n d o m process with a definite scaling factor or a fractal dimension and obeys the inverse power law. Recently West (21) discovered that these systems are much more error-tolerant than the regular nonfractal system, therefore the variability in physiological systems can be much more closely modeled by the fractal system. West and Goldberger were the first to use the fractal concepts to model the bronchial tree (22) emanating from the trachea of mammalian lung, and also the His-Purkinje conduction system of the heart (9) and their results show clearly that these systems have a definite property o f self-similarity which follows very closely to the inverse power law. These types of processes or fractal systems are generally represented in the frequency domain by:
Acknowledgment-The authors wish to express their thanks to the office of Naval Research for its generous support and to Dr. Banu Onaral and Dr. Chris Rorres for their valuable suggestions. Address correspondenceto Dr. H. H. Sun, Department of Electrical and ComputerEngineering, Drexel University, Philadelphia, PA 19104. 597
598
1-1.1t. Sun and A. Charef 1 X(S)
-- -
S- - m
(1)
where s is the complex frequency and m is a positive fractional number and can also be called as the fractal dimension. However, in most cases, the system usually exhibits a finite magnitude at very low frequency or at s--, 0. Hence, we propose that the fractal system which has been represented as the inverse power law equation or 1/f m, where m is any fractional number, can best be represented by a fractional power pole (FPP) as follows:
X(s)
=
1
(2)
where P r has been defined as the corner frequency. This type of expression gives a much wider representation to the natural and physiological phenomena because its low-frequency magnitude is finite instead o f infinite. In addition, the distribution function or the distribution of relaxation time for this type of fractal system also obeys the inverse power law or can be represented by CT ~, where z is the relaxation time. We can show that the time domain representation o f the fractal system as represented by F P P expression is a type of Weierstrass function, a fractal system as originally pointed out by Mandelbrot. In a separate (4) paper we have also developed a singularity method to represent the FPP by a series of pole-zero pairs at the negative real axis, thus the system can be modeled by a series of RC circuits. For more complicated systems where the fractal dimension may be varied from one value to another or systems with multiple fractal dimension, we can then easily represent these systems by a more general equation of multiple FPP, or in this case a higher order fractal system is proposed. To solve the behavior of the system and to predict its output response with different types of input to the system under different initial conditions, we have developed the time domain expression o f the fractal system by a set of linear differential equations with time-varying coefficients. We have made a complete investigation of the time domain representation o f the fractal system in terms o f existence, uniqueness, stability, boundness, and asymptotic behaviour of its solution. These equations will automatically lead to the regular time-invariant system or the nonfractal representation, if the fractal dimensions mi approach to unity or any other integer numbers. We have obtained the exact solution in single fractal system case and also the asymptotic solution of different varieties of multiple fractal systems. These solutions can be easily plotted in terms of impulse response and step response to demonstrate the dynamic behavior of the system under various initial conditions. We then consider the dynamic behavior of systems with a pair of dominant poles whose performance criteria are well-known and we show that the system becomes less oscillatory and also more stable in the fractal domain. Additional examples are presented using the His-Purkinje system for the high order single fractal system and the impedance polarization of inert metal for the multiple fractal system to study their behavior in the time domain.
Fractal System
599 SINGLE FRACTAL SYSTEM-SINGLE ORDER
We shall define the single fractal system or an FPP, as in Eq. 2, where 0 < m < 1 represents the single order case with s = jr and PT is the corner frequency in the log-log Bode plot. The magnitude of X ( s ) at P r is - 3 m dB and the asymptotic slope is zero dB for o~ < PT and - 2 0 m dB/decade for o~ > PT as shown in Fig. 1. In Fig. 1 we have also shown the pole-zero representation (or the singularity function) where Pi and zi are the pole-zero pair in the negative real axis which were obtained by using
PT
vo z0
P1
z~
P2
z~
........
)
Frequency
••-----•
-2Ore.dB/dec]
~~
FIGURE 1A. Bode plot of 1/(1 + S/Pr) m with slope of - - 2 0 m dB/dec and its approximation as zigzag straight lines with individual slopes of - 2 0 dB/dec and 0 riB~dec.
Pc
P0
Zo
P
FIGURE lB. Way to choose the singularities (pole-zero pair) by assuming a prescribed error of y dB between the - 2 0 m dB/dec line and the lines of alternate slopes of - 2 0 dB/dec and 0 dB/dec.
600
H.H. Sun and A. Charef
alternate slopes o f - 2 0 dB/decade and 0 dB/decade with a prescribed error of y dB. Using this approach we have been able to express the F P P as: N
x ( s ) = Y', /=0 (1 + sri)
(3)
where 1
ri
=
(pr)ip ~
Po = the first pole Zi 1)=--
r =
Pi Pi+l - Zi
rp
Pi+l for i = 0,1 . . . . . N. Pi The values of N and the first pole Po are determined by a prescribed error of y dB and the system bandwidth. The details were presented in a separate paper (4). Obviously, i f y ~ 0, then N ~ oo, therefore we can express Eq. 3 as: =
) dr X ( s ) = fo ~ (1G+( r sr)
(4)
where G (r) is the distribution function or the distribution of relaxation time and it is related to Eq. 3 by: N
G ( r ) = ~_j ki6('r - ri) ,
(5)
i=0
where 6(.) is the delta function. This concept was first introduced by Von Schweidler (11) and later Fuoss and Kirkwood (8) derived a formula to obtain the distribution function directly from the original transfer function. Cole-Cole (5) applied this method to find the distribution function for the dispersion and absorption phenomena of a number of liquid dielectrics which were represented as a depressed semi-circle in the Cole-Cole diagram. Later Davidson and Cole (6) found that the complex dielectric for certain types of liquid can be represented by the F P P equation and if m = 1, it becomes the Debye type (7) or a single negative real pole or a single exponential. They also used the FuossKirkwood formula to obtain the distribution function as follows:
G ( r ) = Sin(mTr) ( z r - ~ z )
T~
TT
(6) =0
r>ZT
Fractal System
601
or
G(r) = Co Tm for r << rr ,
(7)
which is true if o~ << ~o, since o~ = 1/r and it is clearly an inverse power law model. I f we take the inverse Laplace t r a n s f o r m o f Eq. 4 o n both sides, we have:
x(t) = s
G(r)r e-'/~ dr
(8)
which is the same as the 1 / f case as described by West and Shlesinger (23) if we replace G ( r ) / r by P(r), where P ( r ) has been defined by them as the distribution function in the non-deterministic case and has been s h o w n that it behaves as the inverse power law. I f we let wi = (k)ico0 where i = 0,1,2 . . . . and k is a constant which can be defined as the ratio o f two successive poles, then ri = r 0 / ( ) Q i and G(ri) is given by: el
G (ri) --- - Or i
(9)
'
where C1 = Corff' and o~ = k m, therefore, G(r) can be written as:
G(T) = k G(Ti)~(T -- Ti) ~- k i=0
i=0
-"7 c l ~(T-- Ti) 9
(lO)
O/t
Substituting G(r) into Eq. 3, which becomes: 1
Ogi
X ( s ) = C1
the inverse Laplace t r a n s f o r m o f this is:
x ( t ) : C 1 ~ ~1 (~ki6oo)e_t~i~o
(11)
i = o o~ t
which is a type o f Weierstrass function with: log o~ m--
log )~
This shows clearly that the F P P system is a fractal system which obeys the inverse power law. Sun et aL (17) has developed a time-varying differential equation to represent the single fractal system as follows:
602
11.1-1. Sun and A. Charef
x+
a+
x=O
,
(12)
where a = P r and b = (1 - m) and it has a solution of: x(t) = ct(m-l)e-tPr .
(13)
E q u a t i o n 11 should therefore converge to Eq. 13 as shown in A p p e n d i x I. We have thus been able to show that the F P P representation o f the fractal system is valid and its time d o m a i n representation is a time-varying differential equation which becomes time-invariant if m = 1.
MULTIPLE FRACTAL SYSTEM, TIME-VARYING LINEAR DIFFERENTIAL EQUATION REPRESENTATION F o r high order single fractal system, we can define the fractal dimension as/3 instead o f rn in Eq. 2 and (N-
1) < / 3 < N ,
(14)
where N is an integer but larger than one. This can be classified as a high order single fractal system. Let us first develop a m o r e general type o f multiple fractal system which can be expressed as follows:
1 X(S)
--
,
0 < m i < 1;
for i = 1,2 . . . . . n
(15)
ffI (s + pi) mi i=1
and this becomes a h i g h order single fractal system if all poles Pi for i = 1,2 . . . . , n are the same and ~ m~ =/3, where/3 is defined in Eq. 14. i=1
Let us define the characteristic equation o f Eq. 15 as the one when all m~ = 1 for i = 1,2 . . . . . n or:
f
i=l
f
-
•
(S + p i ) = sn + al S n - I + a2 Sn-2 +
...
+an
9
(16)
We shall assume for the time being that all the eigenvalues are distinct and the coefficients ak can be easily expressed in terms o f Pi as: ao= 1
ak = -~, i l = l i2=1
"'" Z P i l P i 2 . . .Pik ik=l
(17)
for k = 1 , 2 , . . . , n ; and il ~ i 2 # : 9 9 9 4: ik. W e shall also define r as the coefficients o f the p r o d u c t o f (n - 1) poles for each mi, when i *: k or:
Fractal System
603
ff-~ (s + pk)
mi
= mi(ci,1 Sn-I d- Ci,2 S n - 2 -I- . . .
(18)
-1- Ci, n)
k=l
kr
where Ci, 1 = 1
Ci'k
(k
- 1) v 9
i1=1 i2=1
.
.
.
.
.
i(k_l)=l
for k a n d i = 1,2 . . . . . n; a n d i r i 1 r i 2 r . . . r S i m i l a r to the single f r a c t a l system o f Eq. 12, we can t h u s r e p r e s e n t the m u l t i p l e f r a c t a l system o f Eq. 15 as a n n t h o r d e r linear t i m e - v a r y i n g d i f f e r e n t i a l e q u a t i o n :
(20)
w h i c h is d e r i v e d as Eq. 11.7 in A p p e n d i x II. This can also be r e p r e s e n t e d b y the state v a r i a b l e e q u a t i o n o f the p h a s e v a r i a b l e f o r m as: s = A(t)x + Bu , w h e r e x is the state v a r i a b l e vector, B = [ 0
A(t) =
0
(21) ...
1] r a n d
0
1
...
0
0
0
0
...
0
1
(22) w h e r e ag are d e f i n e d in E q . 17 a n d bk are: n
bk= (n-k+
1)ak_ 1 --
~-~miCik
fork=
1,2 . . . . . . . .
n ,
(23)
i=1
a n d ci, k's are as d e f i n e d b y Eq. 19. T h e detail o f the d e r i v a t i o n is s h o w n in A p p e n dix II. It is interesting to see t h a t Eq. 22 a p p r o a c h e s t o the case o f t i m e - i n v a r i a n t case i f t --, oo o r it c a n be easily s h o w n t h a t if mk = 1 t h e n all b~ = 0, f o r k = 1,2 . . . . n. We c a n also express A ( t ) as: A ( t ) = A1 + A 2 ( t )
,
(24)
604
H . H . Sun a n d A . C h a r e f
where A 1 is a constant n x n matrix which represents the case of mk = 1, for k = 1,2 . . . . . n or the nonfractal case, and A2(t) is a time-varying n x n matrix which becomes the zero matrix as t becomes very large. We can apply this method to the case of a higher order single fractal system as follows:
1 X ( S ) = ( Sn + a l S n _ l + . . .
+ a n - i S + an) m '
(25)
where 0 < m < 1. This represents an nth order single fractal system whose A(t) is defined in Eq. 22 with: b k = (n + 1 -- k)(1 - m ) a k _ l
,
(26)
and a0 = 1 for k = 1,2 . . . . ,n. Similarly, we can also apply this method to a higher order single fractal case with /3 defined by Eq. 14, and we can rewrite the equation as:
1 X(s)
-
[(s + p T ) N ] m '
(27)
where N m = f3 and 0 < m < 1, and N is the lowest integer greater than/5. The characteristic equation becomes an N t h order polynomial in this case. T I M E S O L U T I O N T O T H E M U L T I P L E F R A C T A L SYSTEM The solution of this type of time-varying differential equations have been discussed in detail by Bellman (2). Before we present our solution we shall present a number of theorems on the existence, uniqueness, boundness and asymptotic stability of the solution of this type of time-varying differential equation. The following theorems due to Bellman are presented here without proof.
Theorem 1: Existence and Uniqueness of the Solution Let an n x n matrix H ( t ) be the continuous function in the interval [to > 0, tl ], then in this interval there exists a unique solution of the linear differential equation of
z=H(t)z
(28)
where z(t) is n-dimensional vector, and z(t0) -- z0 is the initial condition vector. F r o m Eq. 22, the n x n matrix A(t) is a continuous function in the interval [to > 0,oo]. Therefore according to theorem 1, given the n-dimensional initial condition x(t0) - xo, there exists a unique solution for Eq. 21 in the interval [to > 0, 0o].
Theorem 2: Boundness Consider the following linear differential equation: i = [H1 + H2(t)]z ,
(29)
Fractal System
605
where z is an n-dimensional vector, n x n matrix. Let
I't I
is an nxn constant matrix, and H2 (t) is an
a) the eigenvalues of H 1 have distinct nonpositive real parts, ~ dH2 b) H2 (t) ~ 0 at t --, 0% and ~ dt < oo
f,
to
where II. II represents the norm of the matrix and under the two above conditions the eigenvalues of [H1 + H2(t)] have distinct non-positive real parts for t >__t~ and t~ is large enough, then all solutions of Eq. 29 are bounded as t ~ ~ . From Eq. 24, it follows that: a) if we assume that all the eigenvalues of the nxn matrix A~ have distinct nonpositive real parts, and b) the matrix A2(t) satisfy the following conditions: A 2 (t) ~ 0 at t ~ 0% and
f
oo dA2 I ---5-7.
to
s
dt < oo
then from the above theorem, we conclude that all solutions of Eq. 21 are bounded as t ~ . Theorem 3: Asymptotic Stability Consider the linear time-varying differential equation of Eq. 29, if all solutions of : nlz
(30)
approaches zero as t ~ 0% then the same is for the solutions of Eq. 29 provided that IIn2(t) II -< C~ for t _> tl, where t~ is large enough, and C1 is a constant which may depend on H1. This proves the asymptotic stability of the solution of Eq. 29 and also of the multiple fractal system. Theorem 4: Asymptotic Form of the Solution Consider Eq. 29 and let: a) the eigenvalues of the matrix H1 to be simple ft ~ dH2 b) H 2 ( t ) - - , 0 a t t ~ , a n d --~- dt
and under the above conditions the eigenvalues of [H1 + H2(t)] are simple for t >_ tl and tl is large enough, then there exist n independent solutions Zk(t) whose asymptotic forms are: t
Zk(t ) ----C k
exp
ft
0
k e ( r ) dr ,
(31)
606
H.H. Sun and A. Charef
where Ck is a constant, nonzero vector, Xk(t) are the eigenvalues of [H~ + Hz(t)], and k = 1,2 . . . . . n. We shall then apply the above theorem to the fractal system defined before and state the following theorem:
Theorem 5: Asymptotic Solution for Multiple Fractal System Consider the following equation:
= [Al+Az(t)]x
,
(32)
where the matrices A1 and A2(t) are defined by Eq. 22 a) if we assume that all the eigenvalues of A1 are simple, and b) the matrix Az(t) = Bt -1 (where B is a constant matrix) satisfy the following conditions:
A2 (t) --, 0 at t ~ oo
and
dt fro ~ dA2
dt <
oo
then, the eigenvalues Xk(t) of [A1 + Az(t)], for t _> tl where tl is large enough, are given as: Xk(t) = - P k + ( m ~ - 1)t -1 ,
(33)
where Pk are the eigenvalues o f the matrix A1 or the nonfractal poles of the system and mk are the fractal dimensions corresponding to each Pk. (Details of the p r o o f is shown in Appendix III.) Therefore, we can obtain the n independent solutions xk(t) whose asymptotic forms are:
Xk(t) = Ckexp
X~(r) d r
xk(t) = Cke-'Pkt~m~-~) ,
(34)
where Ck is a nontrivial vector for 1 _< k _< n.
S U M M A R Y OF SOLUTIONS FOR DIFFERENT TYPES OF FRACTAL SYSTEM We have derived the solution of the single order multiple fractal system. Additional solutions can be derived for various other types of fractal systems. We shall summarize the results by classifying them into different types of fractal system as follows:
1. Single Fractal Simple Order System This is defined in Eq. 2 and it has an exact solution as shown in Eq. 13.
Fractal System
607
2. Single Fractal Multiple Order System This is defined in Eq. 27 with the fractal dimension 13defined in Eq. 14. The exact solution can be easily derived by the use of the above theorems as: x ( t ) = [ ~i=~ lkiti-~
(35)
+ k N t ~ - l ] e-pt
3. Multiple Fraetal Simple Order System This defined in Eq. 15 and the asymptotic solution is given in Eq. 34. 4. Multiple Fractal Multiple Order System The system equation is defined as:
X(s)
-
1
,
q
(36)
I I (s + pi) ~i i=1
where (Ni - 1) < (Ji ( Ni and Ni > 1 for i = 1,2 . . . . . q the asymptotic solution is:
x(t) =
(kj, i t]-I + kj, Nit~i-I)e -pit i=1 [_ j = l
1
,
(37)
which can be easily derived from the above theorems and equations.
D O M I N A N T P O L E S W I T H S I N G L E FRACTAL DIMENSION Numerous studies have been carried out mostly on the representation of the ffactal system to various natural and physiological phenomena in the frequency domain and also about how to obtain the fractal dimension for a single-order fractal system Little has been done about the system with multiple fractal dimensions and the high-order system with single or multiple fractal dimensions. Also, most of the studies have been carried out only on the frequency domain, although many have applied the idea of distribution function of relaxation time, but nothing has been done in the time domain analysis either on time representation or on the time solution of the system with various types of inputs. We shall now try to probe the dynamic behavior of the fractal system and to compare the various performance criteria of the linear nonfractal system with the fractat one. We shall employ the dynamic behavior of a pair of dominant poles (10), which have been used extensively to represent most linear systems as the nonfractal system and compare it with the same system with fractal dimension. Consider the following second-order single fractal system as:
1
X(S)
=
(S 2 +
2~'o0.s +
602) m '
(38)
608
14.1-1. Sun and A . Charef
where ~"is the damping constant and co, is the natural frequency of oscillation of the system F r o m the above discussion, this becomes a linear time-varying second-order differential equation in the time domain as:
J~+
2~'~o.+
t
x = 0
t
(39)
or the matrix A(t) can be represented as: 0 A(t) =
_wz.
1
2~-~on(1 - m )
2~'w.
t
2(1 - m)
t
Solution can then be obtained for the step input or impulse input. Two examples are given as follows: 1. o~n = 25, ~"= 0.4 stable but highly oscillatory case. Figure 2 shows the step response of the system with the fractal dimension varies as 1, 0.7, 0.4, 0.2, 0.05. It is interesting to note that the response becomes less oscillatory as the system takes on the fractal dimension, and the smaller the fractal dimension is the tess
x(O
m=l ,~~/ m=0.7
g
~
t
0.80
i.~o
z.,~'
~'.~
;.oD
FIGURE 2. Step responses of the second order dominant-poles system with fractal dimensions.
Fractal System
609
oscillatory it becomes and it can finally reach the stage of overdamped case with further decrease on fractal dimension. . ~0n = 25, ~"= 0.0, this is the case of sustained oscillations and it is usually considered to be critically stable or unstable. However, if the system takes on the fractal dimension and if the fractal dimension becomes smaller, the response can be considered as a damped oscillation with the final steady state value reaches the zero as limit. In other words, the system becomes stabilized f r o m the case when the system was originally considered to be critically stable under the nonfractal case. This can be easily seen in Fig. 3, where the fractal dimension is varied as 1, 0.6, 0.2. THE H I S - P U R K I N J E SYSTEM
The irregular but self-similar pattern of branchings seen in the His-Purkinje system of the heart has been studied by Goldberger and West. Depolarization of the myocardial cell of the left and the right ventricules via the His-Purkinje tree produces the QRS complex on the surface electrocardiogram. The QRS complex was therefore used to detect this mechanism by analyzing a group of 21 healthy adult males. A plot of the Fourier spectrum is shown in Fig. 4 which shows a log-log plot of the power versus the harmonics. This has been defined as the inverse power law behavior. However, one can easily detect that the system reduced to a constant value at very low frequency, therefore we shall use the F P P representation and obtain the following equation:
xO) m= 1 o
m=.6 T;I=. 2
g
r,i o
<5"
9
1. O0
2.00
3. CO
4. CO
~;. O0
FIGURE 3. Step responses o f a second order critically stable case with different fractal dimensions.
610
H.H. Sun and A, Charef
ORS SPECTRUM
I,
4
9~ -
~
*
*
2
E c~
o
i
32
HARMONICS
i
0.2
I
I
0.6
l
i
1.0
1.4
1.8
log(harmonic) FIGURE 4. The power spectral density of normal depolarization (ORS) waveform from 21 healthy men is shown [from Goldberger et al. (5)]. Reproduced from the Biophysical Journal, 1985, VoI. 48, pp. 5 2 5 - 5 2 8 by copyright permission of the Biophysical Society.
X(s)
-
K (s + p ) ~
'
(41)
w h e r e p = 27r • 4 r a d i a n / s e c as o b t a i n e d at t h e h a l f p o w e r p o i n t a n d fl = 2.1. T h i s is a h i g h o r d e r single f r a c t a l s y s t e m w i t h m = j3/3 = 0.7 o r :
X(s)
-
K [(S q- p ) 3 ]m
T h e A ( t ) o f E q . 22 is a 3 x 3 m a t r i x w i t h t h e e l e m e n t s as: al -- 3 p ,
a2 = 3P 2,
f13 = p 3
and b I = 3(1 - m ) ,
b 2 = 6p(1 - m),
b 3 = 3p2(1 - m ) .
T h e i m p u l s e r e s p o n s e a n d step r e s p o n s e o f t h e s y s t e m a r e s h o w n in F i g . 5. T h i s s h o w s t h a t t h e s y s t e m h a s a v e r y s m o o t h r e s p o n s e a n d it r e a c h e s its z e r o
Fractal System
611
x(t)
2
$
r
~
0.08
09
0.24
0.32
0.40
FIGURE 5A. Impulse response of the His-Purkinje system9
x(t)
o ,-i,
,
~
0.20
!
0.40
09W
0'.80
1'.00
FIGURE 5B. Step response of the His-Purkinje system9
,t
612
H.H. Sun and A. Charef
steady state value at each pulse. The time required to reach its half-life based on the impulse response is about 70 ms. At each pumping of the heart the signal spread out gradually to the numerous branches and finally reduces to zero for each pulse. This can be easily seen from the study which was made on normal healthy adult males. If we increase the slope of the frequency spectrum to be larger than 2.1, then the time required for the half-life becomes larger which shows that the signal will take a longer time to propagate into the tree branches. This might be the case for people with heart diseases or with a rather weak heart. P O L A R I Z A T I O N I M P E D A N C E OF E L E C T R O D E For the multiple fractal system we chose the frequency response of the polarization impedance of Pt electrode as shown in Fig. 6 and as originally developed by Onaral and Schwan (13). The Pt electrode is under small signal excitation or at the linear range, with the geometrical area 0.0855 cm 2, in 0.9~ NaC1 electrolyte solution with pH = 6.5 + 0.1, temperature of 25 + I~ It has been shown that the polarization impedance behaved very much as the inverse power law. Sun et al. (15,16) developed a generalized equation to represent the polarization phenomena by a number of FPP. Recently Tsao et al. (19) developed a computer program to simulate the multiple F P P system, and the result is a three-fractal system as: X(s)
~
I
=
i
(S + pl)ml(s + p2)m2(s + p3) m3 '
i It I
I
I II l
I
I lll'l
I II l
I
! II..
10 ~
10 3
- Io 2 &
102
10
t Ill
(>001
(
O'Ol
i Ill
l
I Ill
l
I ILl.
0'I I 10 fFequency, Hz
i
I Ill
100
i
i tL
1000
FIGURE 6. Series polarization resistance Rp, capacitance Cp, and impedance Zp against frequency. Pt electrode; geometrical areas: 0 . 0 8 5 5 cm 2, linear measurement (65 mV); 0 . 9 % NaCh pH: 6.5 _+O. 1 ; T = 25 _+ 1 ~ [Onaral e t al. (16)]. Reproduced from Medical & Biological Engineering & Computing, 1982, Vol, 20, pp. 2 9 9 - 3 0 6 by copyright permission of the International Federation for Medical and Biological Engineering.
Fractal System
where ml m2 m3 PI P2 P3
= = = = = =
613
0.11 0.36 0.35 0.0001 0.00472 8.065,
and the A(t) of Eq. 22 is 3 x 3 matrix with the elements as: al = P l + P 2 + P 3 a2 = P l P 2 + P i P 3 + PzP3 a3 = P l P 2 P 3
bl=3-(ml+m2+m3) b2 = 2al - [ m l ( p 2 + P 3 ) + m2(Px + P 3 ) + m 3 ( p l +P2)] b3 = a2 -
[ m l p 2 P 3 -I- m z p l P 3 + m 3 p l P z ] .
The impulse-response and the step-response are shown in Fig. 7. It again shows a system with very smooth response and it approaches to its steady state with no oscillations or damping. It took an unusually long time to reach its steady state value and its half-life based on the impulse response is around five minutes. This might be the case for the polarization effect of electrodes and the smooth response is naturally a phenomenon which commonly happens at the fractal system. CONCLUSION The FPP was first introduced by Davidson and Cole (6) to represent the dielectric relaxation phenomenon in Glycerol and Propylene Glycol and to show the difference f r o m the Debye model for n-Propylene. The slopes of the log-log frequency plot shows a range of 0.550-0.666, and the low-frequency value approach to a constant. They also derived the distribution function of relaxation time and show a smooth curve that varies as the inverse power law. Sun et aL (14) have made a study of the F P P and approximates it by using Pade approximation. They also introduced the multiple FPP model to represent the polarization phenomena of noble electrode. Tsao et al. (18) proposed the singularity function approach and designed the R-C circuit model for the system. Recently Charef et al. (4) improved the singularity function method and extended it to the multiple F P P or multiple fractal system and also enabled one to obtain the distribution function for the multiple fractal system. We therefore propose the use of the F P P to represent the fractal system and also prove that the F P P model obeys the inverse power law. The advantages of the F P P model are therefore summarized as: 1. It gives a finite value at the low frequency and enables it to represent the system with this type o f phenomenon.
614
H.H. Sun and A. Charef
x(t)
Ok 0
8
0'
I
r
1. O0
2. O0
,
,
3.00
4.00
5.00
,
t(10%
FIGURE 7A. Impulse response of the polarization impedance of the electrode.
x(o
oo ,-i.
o
~
i 0.20
i 0.40
I O. gO
J
,,.80
!
, t (10 s)
1. O0
FIGURE 7B. Step response of the polarization impedance of the electrode.
Fractal System
615
2. It can be extended to the general model which we have proposed namely, single fractal-single order, single fractal-high order, multiple fractal-single order, multiple fractal-high order. 3, It gives a direct time domain representation in the f o r m of linear time-varying differential equation, whose solutions were available in the asymptotic form. This enables one to study the system responses under different types of input or initial conditions. 4. The F P P is more or less a deterministic model whereas the 1/f is a nondeterministic model. 5. The F P P model can be approximated by the singularity function or a series of pole-zero functions at negative real axis. Thus, one can use the R-C circuit to model the fractal system within a desired error and frequency band. M a n y p h e n o m e n a as in coding and in biochemical processes have been observed to have the fractal behavior been expressed in the form of FPP. Mandelbrot (3) studied the coding of the word in the brain and expressed the number o f words, N ( t ) , or the word's rank, as a function of t, the relative cost, and found that the probability, PN, is directly related to the rank by:
PN--
P ( N + B) ~
where P,B are constants and 1 < 3' < 1.2 for most languages, although in an exceptional case, 3, m a y reach 1.6, for example, in children's talk. The log-log plot o f p N versus N also gives a curve of fractional slope which is an F P P type with constant value at very low N. Similarly, Austin et al. (1) have found that the process of Myoglobin rebinding of carbon monoxide and dioxygen after photodissociation is not an exponential function in time but rather in the fractional power behavior, and can also be expressed in the f o r m of FPP. We have used t h e time domain analysis to obtain the dynamic behavior of fractal system and have found that they are much more stable, less oscillatory and usually have very smooth changes of magnitude in the time domain. West (21) has discovered the error-tolerant nature of the fractal system and we have found this is also true in our R-C circuit model of the singularity representation o f the fractal system. We found that the parameter variation of the R-C circuit has a much greater tolerance and the effect of parameter variation is less sensitive to the pole-zero location. The tree branching nature of the fractal system with its self-similarity property has made the system more error-tolerant and perhaps it also made the system more stable with a better performance than the regular nonfractal system. This is true in most of the biological system where nature has made it more adaptable to the environment and it is more stable and error-tolerant and with smooth and gradual changes instead of abrupt interruption. This is the nature o f the biological system under normal and healthy condition and if these conditions are interrupted, then we have the disturbance and sickness which will m a k e the system no longer adaptable to the environment.
REFERENCES 1. Austin, R.H.; Beeson, K.W.; Eisenstein, L.; Frauenfelder, H.i Gunsalus, I.C. Dynamicsof ligand binding to myhoglobin. 14 (24): 1975.
616
H.H. Sun and A. Charef
2. Bellman, R. Stability theory of differential equation. NY: McGraw-Hill; 1953. 3. Brillouin, L. Science and information theory. 2nd ed. NY: Academic Press; 1962. 4. Charef, A.; Sun, H.H.; Tsao, Y.Y.; Onaral, B. Fractal system, as represented by singularity function. Submitted to IEEE Trans. Auto. Cont. 5. Cole, K.S.; Cole, R.H. Dispersion and absorption in dielectrics, alternating current characterization. J. Chem. Physics 9:341-351; 1941. 6. Davidson, D.W.; Cole, R.H. Dielectric relaxation in glycerol, propylene glycol and n-propanol. J. Chem. Physics 19:1484-1490; 1951. 7. Debye, P. Polar molecules. NY: Dover Pub. 1929. 8. Fuoss, R.M.; Kirkwood, J.K. Electrical properties of solids VIII-Dipole moments in polyvinyle chloride diphenyl systems. J. Ame. Chem. Soc. 63:385-394; 1941. 9. Goldberger, A.L.; Bhargava, V.; West, B.J.; Mandell, A.J. On a mechanism of cardiac electrical stability-the Fractal Hypothesis. Bioph. J. 48:525-528; 1985. 10. Kuo, B.J. Automatic control system 5th ed. Prentice-Hall Inc., 1987. 11. MacDonald, J.R. Impedance spectroscopy. NY: John Wiley; 1987. 12. Mandelbrot, B.B. Thefractal geometry of nature. W.H. Freeman & Co.; 1983. 13. Onaral, B.; Schwan, H.P. Linear and nonlinear properties of platinium electrode polarization, Part 1 frequency dependence at very low frequencies. Med. & Biol. Eng. & Comp. 20:299-306; 1982. 14. Sun, H.H.; Abdelwahab, A.A.; Onaral, B. Linear approximation of transfer function with a pole of fractional power. IEEE Trans. Autom. Con. 29:441-444; 1984. 15. Sun, H.H.; Onaral, B. A unified approach to represent metal electrode polarization. IEEE Trans. Biom. Eng. 30:399-406; 1983. 16. Sun, H.H.; Onaral, B.; Tsao, Y.Y. Application of the positive reality principle to metal electrode linear polarization phenomena. IEEE Trans. Biom. Eng. 31:664-674; 1984. 17. Sun, H.H.; Wang, X.; Onaral, B. Onset of nonlinearity in fractal dimension systems: An application to polarized bioelectrode interfaces. Ann. Biota. Eng. 16:111-121; 1988. 18. Tsao, Y.Y. Fractal concepts in the analysis of dispersion/relaxation process. Ph.D. Dissertation. Drexel University Press; 1987. 19. Tsao, Y.Y.; Onaral, B.; Sun, H. An algorithm for determining global parameters of minimum phase systems with fractal power spectrum. IEEE Trans. Ins. Meas. 38:723-729; 1989. 20. Van Der Zeil, A. On the noise spectra of semiconductor noise and of flicker effects. Physica. 16:359372; 1950. 21. West, B.J. Physiology in fractal dimensions: Error tolerance. Annals of Biom. 1990. 22. West, B.J.; Bhargava, V.; Goldberger, A.L. Beyond the principle of similitude: Renormalization in the bronchial tree. J. Appl. Physiol. 60:1089-1097; 1986. 23. West, B.J.; Shlesinger, M.F. On the inevitability of 1/f noise. Int. J. Modern Physics B3:795; 1989.
I: D E R I V A T I O N O F T H E W E I E R S T R A S S FUNCTION AS A SOLUTION FOR THE FIRST ORDER FRACTAL SYSTEM
APPENDIX
F r o m E q . 11, x ( t ) is g i v e n as:
x ( t ) = C1
~-~ 1
.
--: (~ki~oo)e-t>'~~
i=00L t
t h u s , as )x ~ 1 we s h o u l d h a v e t h e a b o v e e q u a t i o n c o n v e r g e to: x ( t ) = KotCm-1)e -t~r. A s )x ~ 1, we c a n write X = (1 + e ) w h e r e ~ ~ 0, t h e n we c a n m a k e the f o l l o w i n g approximation: Xi = (1 + c ) i =
1 +iE
(I.1)
Fractal System
617
t h e n x ( t ) is given as: x ( t ) = C2
a
e-t~~176
=
C2e-t~~
i=1
~1 e-t~~
(1.2)
i=0 Ot
w h e r e C2 = r
and a = -. k
L e t us define:
u ( t ) = x ( t ) e t~~176 = C2
a
e-t~~176
(I.3)
i=0
f r o m Eq. 1.2, we can write:
x ( k t ) = C1
--=1 (kiwo)e_tXz+l~o -- Cee -t~~ i=0 O/t
-
e-t~~
(I.4)
i=0
a n d also
u(kt)
= x(~kt)e xt~~
=
C2
~
e -t"~~
(I.5)
i=0
A s k ~ 1, kWo ~ Wo ~ cot. T h e n , c o m p a r i n g Eq. 1.3 a n d Eq. 1.5, we have: k u ( t ) = - u ( k t ) + C2 9
(I.6)
Ot
I f we restrict o u r a t t e n t i o n to t h e d o m i n a n t b e h a v i o r o f the s o l u t i o n o f the a b o v e e q u a t i o n , t h e n we can c o n s i d e r that:
(I.7)
u ( t ) = _X u ( k t ) OL
w h i c h has the solution: 0.8)
u(t) = At u
where A can be c o n s t a n t a n d m is d e t e r m i n e d b y direct s u b s t i t u t i o n to be:
,ogo /.t - log X -
log(X)
-
,ogo - -
log X
1 = m -
1
T h e r e f o r e , f r o m Eqs. 1.3, 1.8, a n d 1.9, we can write x ( t ) as: x(t)
= u ( t ) e -t~T = A t m - l e -t~r .
.
(I.9)
618
H.H. Sun and A. Charef A P P E N D I X II: DERIVATION OF L I N E A R T I M E - V A R I A N T D I F F E R E N T I A L E Q U A T I O N FOR M U L T I P L E F R A C T A L SYSTEM
Let's n o w consider the nth order fractal model as represented by: 1
X(s) =
,
O < m i < 1;
fori=
1,2 . . . . . n .
f l (S + pi) m i=l The natural logarithm of both sides of the above equation is:
L n [ X ( s ) ] = ~], - m i L n ( s
i=l
+ Pi)
9
By taking the derivative of both sides of the above expression the following result is obtained:
dX X dX _ -S
~ami(ds I i=l \s+pi/
m2p z + . . . . . . m~ + s + s + pl
I
mn ] ds S + pn
n ~ mi ]-[ (S + pj) i=l j=l
dX _
+
X
ds;
for i C j
(II.1)
X;
(I1.2)
f l (S + Pi) i=l
then,
__ ( s + p i )
i=l
ds
-
i=1
mi__ ( s + p j ) j=l
foriCj.
Let:
f l (s + pi) = ~a ak s"-~
i=1
k=0
2 m/fi (s + i=1 j=l
-- 2 2 k=l i=l
and for i #=j ,
Fractal System
619
where:
ao=l
and
a k = ,-7
. i1=1 i2=1
.
ik=l. .
for k = 1,2 . . . . . n,
ci, a = 1 and
il #:
i2
~:
9 9 9
~: i, ;
, i2=1 ~ " " i ( k _~] ( k - 1 1)! [ ~i~=~ , = l Pi~Pi2...Pi,k-j,]
ci, k -
fork, i=l,2
. . . . . n,
i~il~:iz#:
... ~i(k-1) 9
Therefore, we can express Eq. 11.2 as follows:
k~=oak(Sn--k~s) =--[k=~l i~=lmiCi,k(Sn-kx) ] "
(II.3)
F r o m Laplace t r a n s f o r m properties we k n o w that:
Sq
dX -ds
•{tx (q) + qx (q-l) }
(II.4)
and
sqX = s
(q) }
(11.5)
where
X(q) dqx dt q =
and
q = 0,1,2 . . . . .
T h e n Eq. 11.3 becomes:
k=0
k=0
k=l
i=1
By applying the linearity p r o p e r t y o f Laplace t r a n s f o r m to the a b o v e expression, we have:
s
ak(tx (~-k)) + k=O
~ miCi, k(X (n-k))
(n - k ) a k ( x ("-k-1)) k=O
k=l
i=1
= 0
620
H . H . Sun and A . Charef
Now the Laplace transform of 6(t) is: s s
-
= 1 thus d(1)
- 0
ds
therefore Eq. II.6 can be rewritten as: s
t x (~) +
akt +
( n - k + 1)a~_l -
mici, k
x ("-k)
= s
.
k=l
Taking the inverse Laplace transform of both sides of the above equation, we get the following: tx (n) +
akt +
(n -- k + 1)ak-1 -- ~] miCik
k=l
i=1
X (n-k) =
t6(t)
.
J
For t r O, we have \ (n - k + 1)ae_l - ~
miCi, k
i=l
] x ~n-k) = 6 ( t )
.
(II.7)
t
A P P E N D I X III: DERIVATION OF THE EIGENVALUES OF A (t) Let Xk(t) = - [pk + (1 -- m g ) ] t
for k = 1,2 . . . . ,n
be the characteristic roots of the algebraic equation f ( z ) which is given as:
f(z)
= ~-~ ( Z + Xk(t)) = z n 4 - d l ( t ) z
n-~ + d z ( t ) z n - 2 +
...
+d,(t)
,
(III.1)
k=l
where
dj(t)
=al+
b~t a n d d k ( t ) = a k +
~2
Ck,i
b_~t+ i = o ( / ~ Z i )
for k = 2 , . . . , n
,
with ag and bk are as defined in equations (17) and (23), respectively, for k = 1,2 . . . . . n. F o r t > tl (tl is large enough), the following terms will tend to zero,
k-2 Z Ok,i ----*0 i = 0 l(k-i)
FractalSystem
621
thus, the algebraic equation f(z) is:
f(z) = zn + (al-F ~ )zn-l + (a2 + ~-~)Zn-2+ ... + (an + ~-~) .
(II1.2)
Therefore, the characteristic equation of Eq. 24 approach those of the above algebraic equation, i.e., for t > tl (tl is large enough) the eigenvalues of Eq. 24 are:
~k(t)=--[pk+
(ltmk) ]
fork=l,2 ..... n .
