Journal of The Franklin Institute DEVOTED
TO SCIENCE
AND
THE
MECHANIC
ARTS
Volume 283, Numhw 5
May 1967
A Perturbation Method for the Solution of Linear Matrix DiJkrential Equations by
LOUIS A. PIl?ES
Department of Engineering, University of California, The Aerospace Corporation
Los Angeles, and Consultant,
ABSTRACT: The mathematicalformulation of problems in electric circuit theory, structura engineering, and dynamics by the “state-vector” approach leads to a matrix di$eTential equation of the form dB/dx = A(x)S(x) + C+(X)>.This paper presents a useful method for the approximate
solution
of the matrix
di,erential
equation
boundary conditions. The general procedure is illustrated uniform electrical transmission lines of a general type.
I.
subject
to initial
by applying
and two-point
it to the studg of nm-
Introd7.4ction,
The matrix differential equation @ = A(x)X(s) da
+ 4(x)
I
X(0) = &I
(1)
where X(z) is a coIumn matrix or vector of the nth order, A (2) is a square matrix of the nth order and $(x) is a column matrix or vector of the nth order with the initial condition, S( 0) = X0 arises in the mathematica1 study of a diversified group of physical problems. For example, equations of the type of Eq. 1 arise in the study of nonuniform structural problems, dynamic systems with variable elements, time-varying electrical networks, non-uniform transmission lines and in many other problems in physics and engineering (1-5). When the elements of the matrix A (x) are functions of z of a general type, the solution of Eq. 1 cannot be effected by analytical procedures and the investigator must then resort to methods by which approximate solutions of Eq. 1
357
Louis A. Pipes
may be obtained. In this discussion we present certain procedures of the perturbation type which are useful in obtaining approximate solutions of equations of the form Eq. 1. In the physical problems enumerated above, the elements of the matrices A (5) and Q (2) are known functions of x and, it is desired to determine the vector S(z) by solving Eq. 1 subject to the given indiuE vector 80. II,
P(x)
The Transfer Matrix
The solution of Eq. I may be effected by assuming that the vector S(z) may be represented as the sum of two vectors in the form,
(2)
S(s) = SC(X)+ S,(x) with the initial conditions, S,(O) = X0,
and
S,(O)
= 0.
(3)
The vector X,(x) represents the complementary function of the differential Eq. I and the vector S,(z) may be regarded as the particular integral of Eq. 1. If Eq. 2 is now substituted into Eq, 1, we obtain the following:
dS zc
=
A(z)&(s),
(4)
and
ZP=A(~)+%(x) + 4bJ). In order to solve Eq. 4 we introduce the transfer matrix P(z),
S,(x)
= P(x)&,
P(0)
= I,
P(X)
and write,
= matrix of order (n, n)
(6)
with the stipulation that P(0) = I, the nth order unit matrix. If Eq. 6 is substituted into Eq. 4, the result may be written in the form,
(7) Since the elements of the initial vector SOare arbitrary, Eq. 7 can only be satisfied if dP = A(ci~)P(x), dx The square nth order tram&r matrix P(z) condition P(0) = I.
358
P(0)
= I.
thus satisfies Eq. 8 with the initial
_Toumalof The Franklin Institute
If it is possible to solve Eq. 8 for the transfer matrti P(z), by the following procedure: Let
we obtain X,(z)
S,(x) = P(X)&(~)
(9)
where Q(X) is a column matrix to be determined+ If Eq. 9 is now substituted into Eq. 5, the result is,
OF
The first term of the left member of Eq. II vanishes as a consequence of Eq. 8, and we have
P(x)
2 = &cc).
(12)
By now integrating Eq. 12, we obtain
Q(x) = [= P”(u)+(u)
du.
(13)
If Eq. 13 is now substituted into Eq. 9 the result is
x,(z)
= P(x)
/EWt+$44 du, 0
S,(O)
= 0.
(14)
The vector S(z) is now obtained by substituting Eq. 6 and Eq. 14 into Eq. 2, and thus obtaining X(z)
= P(X) [Ho + [=F(u)~~(u) 0
du].
(15)
This is the formal solution of Eq. 1 with X(O) = X0. Note that once the truns@‘ermat& P(z) is obtained, S(z) is found by the use of Eq. 15. Thus it is apparent that the general solution of Eq. 1 for the vector S(z) depends entirely on being able to obtain the lrunsfer mat& P(z) by solving Eq. 8 with the initial condition P(0) = I. III.
Problems
Involving
Two-Point
Boundary
Conditions
In certain problems of practical importance, the initial vector 80 is not known but information regarding X(z) is given at x = 0 and at 17;= L of the form:
MX(0)
Vol. 283,No. 5, May 1967
+ ivS(I.4) = Q
(16) 359
where M and N are square matrices of order n and G is a column matrix or vector of order n. If we place x = L in Eq. 15 we obtain the vector = P(L)[So
X(L)
w =
+ W] where W is the vector,
(17)
JLP-yu)rb(u) ah.
(18)
0
Therefore, if X(z) as given by Eq. 15 is now substituted boundary condition Eq. 16, the result is MSo + NP(L)[So Equation result :
into the two-point
+ W] = G.
(19)
19 may now be solved for the unknown
SO with the following
So = [M + NP(L)]-l[G
- &Y=(L)W].
(20)
With this for SO, Eq. 15 now gives the solution in the range 0 5 2 2 L. IV.
Computation
of the Transfer Matrix
by a Perturbation
Method
As discussed in Section II, in order to obtain the general solution of the matrix differential Eq. 1, first it is necessary to solve the equation, dP = A(x)P(z)
(211
dX
subject to the initial condition P(O) = I to obtain the &an& mat&x P(x). In the usual application to physical problems, Eq. 21 must be solved in the range 0 5 z $ L. An effective perturbation procedure for solving Eq. 21 consists of first obtaining an nth order square matrix of constants Ao, whose elements are the averages of the elements of the matrix A (x) in the range 0 5 x 6 I;, ao that
Having obtained AO by the use of Eq. 22, we now write A (x) as follows: A(s)
= AO + [A(s)
- AD] = AO + a(x)
(23)
where a(z) is the square matrix of the nth order given by a(x) = [A(z)
360
- AoJ.
(24)
Journ;tl of The Franklin Institute
The matrix a(x) is, therefore, a sort of perturbutiolz matrix about the constant matrix Ao. If the elements of the original matrix A (z) are nearly constant or vary slowly in the range 0 $ z 6 L, then the elements of the matrix a(x) are small relative to the corresponding elements of Ao. If Eq. 23 is substituted into Eq. 21, the equation satisfied by the transfer matrix P(X) is,
g
=
P(0) = I.
A0Pb) + ab)P(x),
(25)
The matrix differential Eq. 25 with the initial condition P(0) transformed (2) into the following integral equation:
[I + Is exp 0
P(x) = exp (AOX) The function exp (&x)
is the m&k exp (AOX) =
( -Aou)a(u)P(u)
= I may be
du],
(26)
exponential junction defined by
2
(AG)~/~c!
(27)
k=O
In order to solve the integral Eq. 26 we try a soIution of the form,
P(x)
=
2 Pk(2).
(28)
b=O
If Eq. 28 is now substituted into Eq. 26 we obtain,
2 pk(x> = exp(A&) p
k-0
+
J:exp ( -Aou)a(u)
g IS(U)
du]. (29)
Equation 29 may be satisfied if we let, (30) (31)
P,(x)
= PO(X)
J’Po( -Ub(u)Pl(u) a
Pk+l(u) = PO(X)
du
k = ~‘ P*(-u)a(u)Pk(u) du, cl
(32)
1, 2, 3,
?? ?? ??
(33)
The sequence of matrix functions P,(s) may be shown (6) to form the uniformly convergent series Eq. 28 under suitable restrictions and is, therefore, the solution of the integral Eq. 26. Vol. 283, No. 5, May 1967
361
Louis A. Pipes
FIG. 1. Non-uniform
V. Application
of the Perturbation
to ccNon-Uniform
Transmission
line.
Method Line
As an interesting and important exampIe of the method described in Section IV, consider the uniform transmission line depicted by Fig. I. The folIowing notation is used to describe the line: (3) x = a space coordinate measured along the line
V(z) 1(x) Z(z) Y(X)
= = = =
S(x)
=
the the the the
complex potential across the line conductors at z current flow at the point o of the line impedance per unit length of the line admittance per unit length of the Iine = the state vector of the line
[II 0 A(x)
-Y(z) P(x)
--Z(x)
=[I
I/ = the “A” matrix of the line
0
= the transfer matrix of the line
The state vector X (5) of the non-uniform line satisfies the following differential equation: (3)
2 = A(zi~)X(z).
(34)
The initial state vector X(0) = SO is now the vector whose elements are the potential Y0 and current 10 at the end z = 0 of the line. The state vector of the line X(o) is related to the trunsfer matrk~ of the line P(z) by the equation, S(x) The transfer matrix P(z)
(35)
= P(x)Xo.
of the line, satisfies the equation dP = A(x)P(x)
(36)
dX
362
Journal
of The Franklin
Institute
Equation 36 will now be solved by the method described in Section XV. We first compute the matrix A0 by the equation,
Ao = i /‘A(s) L 0
0
-zJ
[ -Y*
0
d;z: =
I
(37)
where 20 and YO are the merage values of Z(z) and Y(s) over the length L of the line. We now determine the a(z) matrix by the following equation:
a(z)
= CA(z)
x(x)
= Z(a)
-AI,
= [_;CXi
(38)
-yX’]
where - ZO,
&J(z) = Y(x)
by Eq. 26, the integral equation satisfied by P(a)
-
Yo
(39)
is
P(s) = exp L&a) I + /= exp ( -Aou)a(u)P(u)
1
du .
0
(40)
In this special case we have (3), cash (OX) PO(z) = exp (AG)
-&
i - sinh ( f?z)/Ko
where &
sinh (OX) (41)
= cash (8~)
0 = (ZOYO)II2 = a propagation function = (ZO/YO)~‘~= a characteristic impedance.
1 (42) (43)
The matrix PO(Z) is the transfer matrix of a uniform line whose impedance and admittance per unit length are ZO and Yo, the average values of Z(z) and Y(z) over the length of the line. The matrix P,(r) given by Eq. 31 is the first-order correction matrix which must be added to PO(z) to take into account the variability of the line parameters. If the parameters Z(z) and Y Cx) vary slow@ throughout the length of the line, it us usually sufficient to take
P(s) G PO(Z) t- PI(S) as a good approximation to the true transfer matrix P(s) the first-order correction matrix is P,(x)
= PO(X) I’&(
-u)a(u)Po(u)
(44) of the line. By Eq. 31,
du.
(45)
0
The matrix u(u) is obtained by replacing the letter 2 by u in Eq. 38 and is,
Vot.283.No.5,May 1967
363
Louis A. Pipes therefore, u(u)
-a’j@
= c_YOu,
(46)
I’,( -u) is obtained by replacing x by --u in Eq. 41. Accordingly, multiplication and some aIgebraic reductions, we have 2 /0
a,( -zc)a(u)Po(u)
B(x)
du = [
L-E”(~)
-KID%
C(~IlIKo
-
by direct
+ C(z)1 (47)
--B(x)
where F(a)
= 5 I= cz(u)/& 0
- &y(u)]
cash (20~) du
(48)
B(X)
= i 1% [X(U)/& 0
-
sinh (28u) du
(49)
du.
(50)
C(x) = f
P$Yzc)]
J’C4u)/& + &y(u)] 0
The first-order correction matrix PI(S) by Eq, 45 is cash (0x) PI(S)
as given
1 -&C~b) B(2) +Cb)l 1
- KO sinh (0~)
= -sinh
of the transfer matrix P(z)
cash (0~)
(0X)/& X
[F(x)
-
-B(x)
C(x>]/~~o
’
(51)
If the two matrices of the right member of Eq. 51 are multiplied, the correction matrix PI(Z) is seen to be (52)
Plbl = [;;;;;; ;;::::I where p,,(z)
= B(x)
cash (&IT) + [C(z) - P(X)] sinh (03) sinh (h) - [F(s) + C(z)] co& (0~) JIG ~12b) = Pt C(X)] cash (0~) - B(x) sinh (0~) }&-I P&> = f V(z) ~~~(2) = [F(X) + C(x)]sinh (k> - B(s) cash (es). 2 )
If the variation of the elements of the matrix A(z) is small, then we have I 4~)
364
I <<
20
and
I y(x)
I <<
20
(53) (54) (55) (56)
in the range 0 5 5 S L
for 0 5 x 5 L.
Jaurnal of The Franklin
(57)
Institute
Linear Matrix Di$erential Equations
In such a case the transfer makix P(z) a sufficient degree of accuracy.
is given by P(z)
= PO(X) + PI(S) to
VI. Use of the Transfer Parameter Next we discuss an alternative method which is useful in the solution of equations of the type: (58)
Equation 58 describes the distribution of potential and current in a non-uniform electrical transmission line (3). We now define B(z), the variable propagation function of the line by the equation,
B(x) = [Z(x) Y (x)-p.
(59)
The transfer parameter of the line #(CC) is defined by
95(x) =
/‘S(u) ch.
(60)
D
If Eq. SO is differentiated with respect to X, we obtain
d+ - = B(x). da
Using the relation
d & -d = _*_ - e(x) dx
d$ dx
5 &
(62)
we write Eq. 58 in the form,
-Zdr)] [;I or 0 -Y(x) We now introduce the notation:
-Z(X)/O(X) /O(x)
0
I1v
I[
(63)
. (64)
(65) (66) K ($) is the variable character&k
Vol.283, No. 5, May
1967
impedance of the non-uniform transmission
365
Louis A. Pipes
line. With this notation Eq, 64 may be written as (67) If we Iet and
X=
0
-K
i -l/R
0
B-
then Eq. 67 takes the compact form
1
(68)
(69)
-b=BS.
&
In order to simplify the form of Eq. 69, we introduce the following Iinear transformation :
If we now let
[I
j-p/Z 1 0
V
2:=
(70)
and
M=
i
[
0
K--l/a
(71)
the linear transformation Eq. 70 may be expressed in the form
s =
Alz.
(72)
Now by differentiating Eq. 72 with respect to the transfer parameter 4, and indicating differentiation with respect to 4 by a prime, we obtain S’ =
M'z +
Md.
(73)
Equation 69 may now be written in the form (74)
,.S = ES = BMz. If X’ as given by Eq. 73 is now substituted into Eq. 74, there results M’z + Mz’ = BMz.
(75)
We now solve Eq. 75 for .z’, and obtain 2’ =
366
[i,k13M
-
(76)
JP&-‘-J.
Joumal of The
Franklin
Institute
0
1t can easily be shown that
-1
= [ -1
M-IBM
and
0
1
(77)
(78)
If we now substitute Eqs. 77 and 78 into Eq. 76, we have -K’/2K
-1
er
In many cases of practical importance, it happens that the characteristic impedance K varies slowly throughout the range 0 5 x 6 L. In such cases, we have K’/2K 1 < 1. - (30) If we now define Co and C(9) by the equations,
/I 1 0
co=
-Kf/2K
-1
-1
0
ct+>
0
= 0
K’/2K
1 (811
we may then write Eq. 79 in the form,
$
= coz + C(l$)z.
(82)
Equation 82 is similar to Eq. 25 in form, therefore, it can be transformed into the integral equation, (33) where z(0) is the initial vector z at up= 0. In this case, we have (3)
co& (8) exp (Cd)
-sinh where u
Vol. 283, No. S, May
1967
-sinh
(4)
=
and
(4)
cash (#)
1
(34)
0(u) = [Z(U)Y(U)]~‘~.
367
Louis A. Pipes
The initial vector x(0) may be obtained by inverting Eq. 72 and placing 4 = 0, so that K-l’Z(O) 0 x(0) = .X(O) = M-l(O) S(0). (85) 0 Kl@ (0) i 1 The first term of the right member of integral Eq. 83 is
a(d
= exp (C04)2(0) = exp (GdJ-F(O)X(O).
(86)
Let the state vector S(4) that corresponds to the approximate solution x0($) be denoted by SO(4). We then have, as a consequence of Eqs. 72 and 86,
X0($> = IV($) exp (CO+>M-~(O)X(O). If the indicated matrix multiplications may be written as follows:
in Eq. 87 are performed,
(87) the result
L-[K(z)K(O) J-“z L
[K(z)/K(O)
J cash 9
-[K(x)K(O)]~/~
sinh 4
[K (0) /K(z)
sinh 4
-Jcash 0
(@I
where e(2) = C-z(z) U(z)-J1’2,
K(z)
= Z(z)/a(s).
The result, Eq. 88, is the W.K.B.J.L. (Wentsel-Kramers-Brillouin-JeffriesLiouville) approximation (7) . Note that the W.K.B.J.L. approximation reduces to cash (es) (89) -sinh
(6$)/K
in the simple cam of a uniform line in which the functions Z(X) and Y(z) become constants. The W.K.B.J.L. approximation is obtained by using onIy the first term of the right member of integrd Eq. 83. This integral equation has a solution of the form (90)
a(4) = &l(4) T&=0 where
(91)
.%i(+) = exp (cbQ)40) &I($) = exp (Cd) J’ exp (-Cd~)C(u)x~~(u) 0 368
du
n = 1, 2,3, ‘*a.
(92)
Journal of The Franklin Institute
Linear &la&ix DifferentialEquations The state-vector of the line rS, is given by
x = M(4)z($).
(93)
It can therefore be seen that the term ZO($) in Eq. 90 leads to the W.K.B.J.L. approximation. The terms zn($), n = 1, 2, 3, 4, lead to more accurate results than those obtained by the use of the W.K.B.J.L. method. The vector ZI($) is given by Eq. 92, in the form ?? ??
Q(4)
=
exp (C&) /“” exp (-C~U)C(U)ZO(U) 0
clu.
(94)
In cases where the chtiacteristic impedance K(X) varies slowly throughout the range 0 5 z 5 L, suflkient accuracy is usually obtained by the improved W.K.B.J.L. approximation, given by 8 = M(d)ko(&) VII.
Transmission
Non-Uai$wm
a(4>1.
+
Line with Constant
w Propagation
Function
Many non-uniform transmission lines of importance in practice have the property that their propagation function O(X) is a constant, so that 13(z) = cZ(z) Y(z)]~‘~
= 0
(a constant).
(96)
In these special cases, the lrccnsferparameter 4 discussed in Section VI is given by l#l = Jku 0
= ex.
(97)
In this case, it is convenient to write Eq. 58 in the following form: 0
-Z(z)/0
-Y(z)/0 where K(z)
0
0
2 0 I
/I -l/K(s)
--K(z) 0
I[1(981 Fi
I
is the variable characteristic impedance of the line given by K(x)
=
Ei$ = A-_ Y(x) *
In the notation of Eq. 68, Eq. 98 may be rewritten as
Vol. 283,No. 5, May 1967
369
Louis A. Pipes
of we again introduce the transformation Eq. 71, S = Mz, it is easy to show that &‘q. 100 is transformed into the equation,
By direct matrix multiplication,
0 (Ill-lSM>
= [ “,
we obtain
(ilk1%)
;I,
= i”‘?
__KU12K].
(102)
Therefore, Eq. 101 may be written in the form,
(103)
If we now let Nj”,
;],
C(X) =[-K;
K,;2J.
(104)
With this notation, Eq. 101 may be written as dx - = Nz + C(X)& dx
(105)
The integral equation form of Eq. 105 is, z=
exp (N%I /z(O) + J’
exp
(-Nu)c(u)z(u)
du
We have (5) cash (0X) exp(ivcc)
The solution of Eq. 106 is,
-sinh
(kc)
=
L -sinh
(106)
. I
0
(8~)
cash (0~)
I-
(107)
z(2) = 2 G.(x) 9-o
(108)
with .20(z) = exp (iv~)x(O)
(109)
and
znb) = exp (Nz) J
r* ~XP (-Nu)C(u)z,_~(u)
du,
n = 1, 2, 3, ...
(110)
0
370
~ournsl of The
Franklin
Institute
~~(2)
is given by the equation,
a(x) = exp (Nx) [‘exp 0
(-Nu)C(u)
(111)
exp (Nzc) &.X(O).
If we let P(s)
= K’(z)/2K(s)
exp (-Nu)C(u)
F(zc). 1 (112)
we have, - cash (28~)
sinh (20~)
- sinh (2824)
cash (28~)
(113)
exp (iVU) =
We now let,
(114) (1151 We may then write Eq. 111 in the expanded
form,
-= x(O)
(&)J(z)
- cash (~x)H(z)]
[cash (&)J(z)
[sinh
&)H(x)
- cash (Bz)d(x)]
[H(X) cash (8~) - sinh @XV(X)]
If 1 P(X) j = 1 K’(x)/~K(x) solution of Eq. 106 is
/ << 1, in the range 8(X) = z,(2)
The solution
+ G(X).
- sinh (&)A(x)]
* 1 (116)
[sinh a(x)
0 2 x $ L, a satisfactory (117)
of Eq. 98 is then given by,
where x(0)
= M-LAY(O).
(119)
Bibliography (1) L. A. Zadeh and C. A. Desoer, ‘
The State Vector. Approach,” New York, McGraw-Hill Book Co., 1983. P. M. DeRusso, R. J. Roy, and C. M. Close, “St& Variables for Engineers,” New York, John Wiley and Sons, Inc., 1965. L. A. Pipes, “Direct Computation of Transmission Matrices of Electrical Transmiesioo Lines,” Parts I and II, Jou~. Frank. Inst., Vol. 281, Nos. 4 and 5, April and May 1966. B. I<. Kinariwals, I.R.E. Convention Record, Pt. 4, pp. 268-276, 1961. L. A. Pipes, “Matrix Methods for Engineers,” Englewood Cliffs, N. J., Prentice-Hall, Inc., 1963. F. G. Trieomi, “Integral Equations,” New York, Interscience Publishers, 1857. J. C. Water, “Microwave Transmission,” New York, McGraw-Hill Book Co., 1942.
Vol. 283,No. 5, May 1967
371