- Email: [email protected]
An economical method for the numerical solution of convolution equations
B. S. Pariiskii
208
REFERENCES 1.
POLYAK, 8. I’., Iterative methods using Lagrange multip~ers for solving extremal problems with constr~ts of the equation type. Zh. v2hisf. Mat.mat.Fiz., 10, 5, 1098-1106, 1970.
2.
POLYAK, B. T., On the acceleration of the convergence of the convex programming method. Dokl. Akud. Nuuk SSSR, 212,5,1063-1066,1973.
3.
ROBINSON, S. M., A quadratically-convergent Math. Progr4rn.. 3,2, 145-156, 1972.
4.
FADDEEV, D. K. and FADDEEV, V. N., Computational methods of hear ulgebra (Vychislitel’nye metody lineinoi algebry), Fizmatgiz, Moscow-Leningrad, 1963.
5.
COHEN, A., Rate of convergence for root finding and optimization California, Berkeley, 1970.
6.
SMOLYAK, S. A., A qua~atic convergence of the method of conjugate gradients. Proceedings of the Third WinterSchool on ~4thern4fic~l Programming (Tr. III Zimnei shkoly po matem. programmirovaniyu), 558-574, MISI, Moscow, 1970.
I.
MAISTROVSKII, G. D., Proof of the quadratic convergence of the method of conjugate gradients. In: Computational mathematics and computing techniques (Vj%hisl. Matem. i vychl. tekhn.) II, 3-6, FTINT Acad. Nauk Ukr SSR, Khar’kov, 1971.
AN
algorithm for general non-linear programming problems.
algorithms. Ph. D. Dim., Univ.
ECONOMICAL METHOD FOR THE NUMERICAL SOLUTION OF CONVOLUTION EQUATIONS* B. S. PARIISIUI MOSCOW
(Revised 15 January 1976)
A NUh4ERICAL method of solving Fredholm integral equations of the second kind with a kernel dependent on the difference of arguments is presented. The method requires of the order of n memory cells of a computer and n2 arithmetic operations to obtain a solution with an error of the order of nm4. We consider the reduced Wiener-Hopf equation f cP(G+h
s
~(s-~)~(s)~~=~(~),
0
(0
0
the first two derivatives of the kernel and r&&t side exist and be bounded (K(y) may have a discontinuity at y = 0, as usually happens for K(ly I)). Th en, using the formula for tangents, accurate to quantities of the order of n-2 this equation can be replaced by a system of linear algebraic equations with the Toeplitz matrix T@=F, or Let
(2) where
*Zh. @hisI. Mat. mat. Fis., 17, 2,501~504,
1977.
209
Short communications
The resulting system (2) can be solved by well-known economical methods requiring of the order of 4n cells of the memory of a computer and 4n2 operations (multiplications and additions) for the symmetrical case and rather more for the general case [l-5] . Here the methods of [3-51 require modification; it is not necessary to find the inverse matrix T-l as the authors propose; it is only necessary to find triangular matrices P, Q such that PTQ=E, and then find the solution of the system. Here there is no need to store the whole matrices P, Q, obtained row by row or column: since @=PIF=QPF=QY, Y =PF, then finding the j-th row of P, we can find +j =
Pjifi, ,=”
(the j-th component of the vector \k) and then knowing the j-th column of Q (in the symmetrical case it is identical with the j-th column of P, apart from a constant multiplier) in all the components pk of the vector @ we find the next,j-th term of the sum n (pk =
qkj$j. c J=k
In this modification the method described in [3, S] is more economical in computing time than all the others. We present this algorithm:
r=(l+k,,)-I,
ct=rki,
cll=l, i=l,
cpi=0,
(Po=fo,
&=Eo,,=~oo=I,
i =
2:. . . , n;
+2 (...,--
*I,
/So=::,,
a0=ci,
~_~j=~_lj=~_ij=O,
;
j=O, 1, . . . n;
forj=1,2,...,nweput &,j-l=f;i,j-i,
i=o,
j-l,
I,...,
i=j, j-1,. . . ,O;
~ij=Ei-i.j--i-Uj-l~j-l-i,j-1, qtj=rli-l
i=i,
j-l-_j-It;,-l-i,j-1,
aj =
j-l,.
. . ,o;
(i EijCi+l)
Aj-I,
i=O
pj=
(i’liiC-ii+*i)
Qj= (iqijfi
Aj-‘,
i=O
(Pi”(Pi+EijQj,
) Aj_i;
i=o
i=o,
1,. . .
)
j;
Finally ‘Pi=‘q)i,
i=O, 1,. . . , n.
The algorithm is obviously simplified in the symmetrical case. We note that dispensing with the array l, has by [4] increased the computing time (ALGOL, BESM-6). In the methods described in [ 1,2] and the modified methods of [3-5] the rounding errors in the control calculations were negligibly small right up to n = 1000. The necessary and sufficient condition for the application of these methods is the solvability of Eq. (1) with an upper limit of integration z, Osz~l. This is satisfied, for example, if the operator on the left side of the equation is positive-definite, or h= a+iP, p+o, (K(j) is a real function (in the general case all the quantities may be complex). The convergence of the method and its stability to errors in the initial data follows from the results of [6], obtained by using the rectangle formula (our case
B. S. Pariiskii
210
differs only by a shift of the right side, hence all the estimates of [6] hold). Denoting the original operator by L, the approximate one by Lh and the projection operator by Ph ( PhU=U at the points h/2, 31t/2,. . . , (n+l/2)h) ), we have (3) provided that bounded derivatives f and f’, and also K’ and K” exist (KCy) may have a kink at y = 0). For smootherfand K the solution qh can be improved by Runge’s rule, if we also have the solution q3h. @=f
Indeed, let us consider the equation it. We estimate the quantity Ah=Lh (PhCp)
and the equation k[h=Phf=fh In our case we have
approximating
7% I LqJ=cp(t) +h K(s-t)cp(s)ds, J i=o -P&p.
Lh(ph-(ph
ki-jqhi.
j -b h
E
Rutting g(t, s)=M(t-s)cp(s),
wi obtain
n
A,=+’ where
c
g” +
+h2&,(g’(f, t+o) --g’(t,
t-o))
+O(h4)=Phvh2+0(h’),
0
1 u(t) =
G
’
J
g”
as+
+@.
t+o1 -g’(t, t-o)),
0
g”(t, t) = -+t,
t+O) +g”(t, GO)).
We have assumed that bounded derivatives of g up to the fourth order exist for t # s, for which the existence of such derivatives of the kernel and right side is sufficient (we again assume that the kernel has a kink!). We then have Lh(Ph(F-(ph)
=
(Lhph(p-P&)
+
(Ph@-Lh(ph)
=&+Phj-jh=PhVh2+O(h4).
On our assumptions bounded derivatives u’, Y” exist, therefore taking into account (3), where f has been replaced by V,we obtain Phcp-cF,I=h*PhL-‘v+O(h~). Correspondingly, for the step H = 3h we have From this ((ph-(P,~)/8=h’PhL-‘L’+O(h*),
Phq=(ph+
((Ph-(PII)/8+0(h4)
at those points where VH is defined (cphis also defined there, and
Pfl=Ph(p).
A calculation with doubled n makes it possible to obtain similarly (at the centres of the segments of the scale) the improvement (ph.h/z of the linearly interpolated solution &, , and also to improve the derivative and integral of @& ; in all cases we obtain an accuracy of order h4 (for a sufficiently smooth kernel and right side of the equation); &.h/z for not very large third and higher derivatives of cpis more accurate than cpn,h/s, obtained by the methods indicated above from oh and (phi3-
211
Short communications
A more detailed treatment of topics similar to those considered above, the corresponding programs in ALGOL and some new economical methods of solving such systems and integral equations can be found in [7]. ‘Ihe author thanks A. A. Abramov for useful advice. Translated by J. Berry. REFERENCES 1.
CHEPURINA, I. V., Solution of systems with a Toeplitz matrix in the packet of the linear algebraSA3lA. In: ~ume~cal analysis by FORTmN (Chislennyi analiz na FORTRANe), No. 3,152-162, Izd-vo, MGU, Moscow, 1973.
2.
TREITEL, S. and ROBINSON, E. A., High-resolution digital filters. lfX5 Trans. Geosc. Electronics, 4, 1,36-38.1966.
3.
TRENCH, W. F., An algorithm for the inversion of finite Toeplitz matrices. Industr. Appt, 12,3,515-522,1964.
4.
KUTIKOV, L. M., The structure of matrices which are the inverse of the correlation matrices of random vector processes. Zh vi;chisl. Mat. mat. Fiz., 7,4, 764-773, 1967.
5.
JUSTICE, J. H., An algorithm for inverting positive definite Toeplitz matrices. J. AppL Moth., 23, 3, 289-291,1972.
6.
OSTROWSKI, Sur ~approx~ation du d~te~inant de Fredholm par les d~te~in~ts d’equations lineaires.Arkiv Mat., Astmn. och Fys., 26,4, l-1.5,1939.
7.
PARIISKII, B. S., Economical methods for the numerical solution of convolution equations and systems of algebraic equations with Toeplitz matrices (Ekonomichnye metody chislennogo resheniya uravnenii v svertkakh i sistem algebraicheskikh uravnenii s teplitsevymi matritsami). Vts Akad. Nauk SSSR, Moscow, 1976.
des systemes
THE PROBLEM OF IDENTIFYING THE FUNCTIONS GENERATING A NON-LINEAR OPERATOR* R. D. BAGLAI Novosibirsk (Received 11 July 1975; revised lONovember 1975) AN INVERSION formula for an integral equation is given and the problem of identifying four functions generating a non-linear Nemytskii operator is solved. In the investigation of the physical properties of the elements of radioelectronics, optics, electrotechnics the problem arises of identifying the functions V, Q, SEC*, generating the non-linear operator
( f m(W)+$ (q(t))=
V(cp(O)+Q
1’(v).
In [l] it was shown that the unique continuously differentiable solution of problem (1) subject to the conditions V(O)=Q (0) =S (0) =0 can be obtained for known values of the operator Ton the two-parameter family of signals {VP= z sin ot}. Formulas defining the unknown functions were also given there. *Zh. vychisl. Mat. mat. Fiz., 17, 2,504-507,
1977
208
REFERENCES 1.
POLYAK, 8. I’., Iterative methods using Lagrange multip~ers for solving extremal problems with constr~ts of the equation type. Zh. v2hisf. Mat.mat.Fiz., 10, 5, 1098-1106, 1970.
2.
POLYAK, B. T., On the acceleration of the convergence of the convex programming method. Dokl. Akud. Nuuk SSSR, 212,5,1063-1066,1973.
3.
ROBINSON, S. M., A quadratically-convergent Math. Progr4rn.. 3,2, 145-156, 1972.
4.
FADDEEV, D. K. and FADDEEV, V. N., Computational methods of hear ulgebra (Vychislitel’nye metody lineinoi algebry), Fizmatgiz, Moscow-Leningrad, 1963.
5.
COHEN, A., Rate of convergence for root finding and optimization California, Berkeley, 1970.
6.
SMOLYAK, S. A., A qua~atic convergence of the method of conjugate gradients. Proceedings of the Third WinterSchool on ~4thern4fic~l Programming (Tr. III Zimnei shkoly po matem. programmirovaniyu), 558-574, MISI, Moscow, 1970.
I.
MAISTROVSKII, G. D., Proof of the quadratic convergence of the method of conjugate gradients. In: Computational mathematics and computing techniques (Vj%hisl. Matem. i vychl. tekhn.) II, 3-6, FTINT Acad. Nauk Ukr SSR, Khar’kov, 1971.
AN
algorithm for general non-linear programming problems.
algorithms. Ph. D. Dim., Univ.
ECONOMICAL METHOD FOR THE NUMERICAL SOLUTION OF CONVOLUTION EQUATIONS* B. S. PARIISIUI MOSCOW
(Revised 15 January 1976)
A NUh4ERICAL method of solving Fredholm integral equations of the second kind with a kernel dependent on the difference of arguments is presented. The method requires of the order of n memory cells of a computer and n2 arithmetic operations to obtain a solution with an error of the order of nm4. We consider the reduced Wiener-Hopf equation f cP(G+h
s
~(s-~)~(s)~~=~(~),
0
(0
0
the first two derivatives of the kernel and r&&t side exist and be bounded (K(y) may have a discontinuity at y = 0, as usually happens for K(ly I)). Th en, using the formula for tangents, accurate to quantities of the order of n-2 this equation can be replaced by a system of linear algebraic equations with the Toeplitz matrix T@=F, or Let
(2) where
*Zh. @hisI. Mat. mat. Fis., 17, 2,501~504,
1977.
209
Short communications
The resulting system (2) can be solved by well-known economical methods requiring of the order of 4n cells of the memory of a computer and 4n2 operations (multiplications and additions) for the symmetrical case and rather more for the general case [l-5] . Here the methods of [3-51 require modification; it is not necessary to find the inverse matrix T-l as the authors propose; it is only necessary to find triangular matrices P, Q such that PTQ=E, and then find the solution of the system. Here there is no need to store the whole matrices P, Q, obtained row by row or column: since @=PIF=QPF=QY, Y =PF, then finding the j-th row of P, we can find +j =
Pjifi, ,=”
(the j-th component of the vector \k) and then knowing the j-th column of Q (in the symmetrical case it is identical with the j-th column of P, apart from a constant multiplier) in all the components pk of the vector @ we find the next,j-th term of the sum n (pk =
qkj$j. c J=k
In this modification the method described in [3, S] is more economical in computing time than all the others. We present this algorithm:
r=(l+k,,)-I,
ct=rki,
cll=l, i=l,
cpi=0,
(Po=fo,
&=Eo,,=~oo=I,
i =
2:. . . , n;
+2 (...,--
*I,
/So=::,,
a0=ci,
~_~j=~_lj=~_ij=O,
;
j=O, 1, . . . n;
forj=1,2,...,nweput &,j-l=f;i,j-i,
i=o,
j-l,
I,...,
i=j, j-1,. . . ,O;
~ij=Ei-i.j--i-Uj-l~j-l-i,j-1, qtj=rli-l
i=i,
j-l-_j-It;,-l-i,j-1,
aj =
j-l,.
. . ,o;
(i EijCi+l)
Aj-I,
i=O
pj=
(i’liiC-ii+*i)
Qj= (iqijfi
Aj-‘,
i=O
(Pi”(Pi+EijQj,
) Aj_i;
i=o
i=o,
1,. . .
)
j;
Finally ‘Pi=‘q)i,
i=O, 1,. . . , n.
The algorithm is obviously simplified in the symmetrical case. We note that dispensing with the array l, has by [4] increased the computing time (ALGOL, BESM-6). In the methods described in [ 1,2] and the modified methods of [3-5] the rounding errors in the control calculations were negligibly small right up to n = 1000. The necessary and sufficient condition for the application of these methods is the solvability of Eq. (1) with an upper limit of integration z, Osz~l. This is satisfied, for example, if the operator on the left side of the equation is positive-definite, or h= a+iP, p+o, (K(j) is a real function (in the general case all the quantities may be complex). The convergence of the method and its stability to errors in the initial data follows from the results of [6], obtained by using the rectangle formula (our case
B. S. Pariiskii
210
differs only by a shift of the right side, hence all the estimates of [6] hold). Denoting the original operator by L, the approximate one by Lh and the projection operator by Ph ( PhU=U at the points h/2, 31t/2,. . . , (n+l/2)h) ), we have (3) provided that bounded derivatives f and f’, and also K’ and K” exist (KCy) may have a kink at y = 0). For smootherfand K the solution qh can be improved by Runge’s rule, if we also have the solution q3h. @=f
Indeed, let us consider the equation it. We estimate the quantity Ah=Lh (PhCp)
and the equation k[h=Phf=fh In our case we have
approximating
7% I LqJ=cp(t) +h K(s-t)cp(s)ds, J i=o -P&p.
Lh(ph-(ph
ki-jqhi.
j -b h
E
Rutting g(t, s)=M(t-s)cp(s),
wi obtain
n
A,=+’ where
c
g” +
+h2&,(g’(f, t+o) --g’(t,
t-o))
+O(h4)=Phvh2+0(h’),
0
1 u(t) =
G
’
J
g”
as+
+@.
t+o1 -g’(t, t-o)),
0
g”(t, t) = -+t,
t+O) +g”(t, GO)).
We have assumed that bounded derivatives of g up to the fourth order exist for t # s, for which the existence of such derivatives of the kernel and right side is sufficient (we again assume that the kernel has a kink!). We then have Lh(Ph(F-(ph)
=
(Lhph(p-P&)
+
(Ph@-Lh(ph)
=&+Phj-jh=PhVh2+O(h4).
On our assumptions bounded derivatives u’, Y” exist, therefore taking into account (3), where f has been replaced by V,we obtain Phcp-cF,I=h*PhL-‘v+O(h~). Correspondingly, for the step H = 3h we have From this ((ph-(P,~)/8=h’PhL-‘L’+O(h*),
Phq=(ph+
((Ph-(PII)/8+0(h4)
at those points where VH is defined (cphis also defined there, and
Pfl=Ph(p).
A calculation with doubled n makes it possible to obtain similarly (at the centres of the segments of the scale) the improvement (ph.h/z of the linearly interpolated solution &, , and also to improve the derivative and integral of @& ; in all cases we obtain an accuracy of order h4 (for a sufficiently smooth kernel and right side of the equation); &.h/z for not very large third and higher derivatives of cpis more accurate than cpn,h/s, obtained by the methods indicated above from oh and (phi3-
211
Short communications
A more detailed treatment of topics similar to those considered above, the corresponding programs in ALGOL and some new economical methods of solving such systems and integral equations can be found in [7]. ‘Ihe author thanks A. A. Abramov for useful advice. Translated by J. Berry. REFERENCES 1.
CHEPURINA, I. V., Solution of systems with a Toeplitz matrix in the packet of the linear algebraSA3lA. In: ~ume~cal analysis by FORTmN (Chislennyi analiz na FORTRANe), No. 3,152-162, Izd-vo, MGU, Moscow, 1973.
2.
TREITEL, S. and ROBINSON, E. A., High-resolution digital filters. lfX5 Trans. Geosc. Electronics, 4, 1,36-38.1966.
3.
TRENCH, W. F., An algorithm for the inversion of finite Toeplitz matrices. Industr. Appt, 12,3,515-522,1964.
4.
KUTIKOV, L. M., The structure of matrices which are the inverse of the correlation matrices of random vector processes. Zh vi;chisl. Mat. mat. Fiz., 7,4, 764-773, 1967.
5.
JUSTICE, J. H., An algorithm for inverting positive definite Toeplitz matrices. J. AppL Moth., 23, 3, 289-291,1972.
6.
OSTROWSKI, Sur ~approx~ation du d~te~inant de Fredholm par les d~te~in~ts d’equations lineaires.Arkiv Mat., Astmn. och Fys., 26,4, l-1.5,1939.
7.
PARIISKII, B. S., Economical methods for the numerical solution of convolution equations and systems of algebraic equations with Toeplitz matrices (Ekonomichnye metody chislennogo resheniya uravnenii v svertkakh i sistem algebraicheskikh uravnenii s teplitsevymi matritsami). Vts Akad. Nauk SSSR, Moscow, 1976.
des systemes
THE PROBLEM OF IDENTIFYING THE FUNCTIONS GENERATING A NON-LINEAR OPERATOR* R. D. BAGLAI Novosibirsk (Received 11 July 1975; revised lONovember 1975) AN INVERSION formula for an integral equation is given and the problem of identifying four functions generating a non-linear Nemytskii operator is solved. In the investigation of the physical properties of the elements of radioelectronics, optics, electrotechnics the problem arises of identifying the functions V, Q, SEC*, generating the non-linear operator
( f m(W)+$ (q(t))=
V(cp(O)+Q
1’(v).
In [l] it was shown that the unique continuously differentiable solution of problem (1) subject to the conditions V(O)=Q (0) =S (0) =0 can be obtained for known values of the operator Ton the two-parameter family of signals {VP= z sin ot}. Formulas defining the unknown functions were also given there. *Zh. vychisl. Mat. mat. Fiz., 17, 2,504-507,
1977