- Email: [email protected]
Integral equations for potential problems with the source function not located on the boundary
Computers & Smcturcs,
Vol. 8. pp. ll?-115.
Pergamon Press 1978.
Printed in Great Britain
INTEGRAL EQUATIONS FOR POTENTIAL PROBLEMS WITH THE SOURCE FUNCTION NOT LOCATED ON THE BOUNDARY G. DE MEY
Laboratory of Electronics, Ghent State University, Sint Pietersnieuwstraat 41, 9000 Ghent, Belgium (Received December 1975)
Abstract-A special integral equation method for solving Laplace equation in a two dimensional area is studied. Usually an unknown source function is defined along the boundary. In the method to be outlined here, this source function is located on another closed curve. This new method offers several advantages as compared to the classical integral equation technique.
1.INTRODUCTION
Many problems in physics and engineering which are usually described by partial differential equations and appropriate boundary conditions can also be described by equivalent integral equations [l-5]. Integral equations have proved to be very useful for the numerical solution of a given problem&161. If the original problem is two dimensional, the corresponding integral equation will be one-dimensional. This fact reduces the complexity of the computer programming and the computation time is also diminished. For the particular case of a two dimensional potential problem described by the Laplace equation: V’4 = 0
(1)
in an area S with boundary C, the potential 4 is written as: 4(f) = f p(i’)G($‘) dC
(2)
where G(ili’)=llnL
2Tr JJ- F”)
is the Green’s function of the eqn (l)[l]. p(fi is an unknown source function defined along the boundary C. By imposing the boundary conditions on the proposed solution (3), one obtains an integral equation for the function p(fl. Note that the Laplace eqn (1) has been replaced by a one-dimensional integral equation. In the expression (2), the function p(F) is defined along the boundary C of S. However, there is no a prior’ reason why the function p(!) should be defined along C. Taking a second closed line C,, which completely surrounds the area S, one can write the potential 4 as:
the unknown function p(Q, which is now defined along C,. If C, does not coincide with C, the fundamental discontinuities for a boundary condition such as V&E” = 0 should not be taken into account[l7]. One can simply write that: V4(n =
f
p(i')G(i/i') dC;.
(4)
Because the Green’s function (3) is a particular solution of (I), the expression (4) is also a solution of Laplace equation. We have now to impose the boundary condition on (4), in order to obtain an integral equation for CAS Vol. 8, No. 1-H
p(Y)VG(@‘) dC;
(5)
where (5) is still valid on the boundary C. Another feature caused by introducing C, is related to the numerical solution of the integral equation. The continuous function p(?) is then replaced by a set of n numbers pI . . . p. along C,. If C, coincides with C, the numerical noise will be large in the vicinity of the boundary, because the numerical discretisation process takes place there. If C, does not coincide with C, the numerical noise can be made sufficiently low provided C, is not too close to C. This problem will be investigated in this paper for a particular example where C1 is a circular and C a rectangular boundary shape. The boundary conditions on C are chosen in such a way that the problem can also be solved analytically. The accuracy of the numerical results can then easily be checked. 2. ExAMF%E
The numerical solution of Laplace equation in an area S with boundary C by means of an integral equation with the source function defined along a closed line C, (Z C) will be given for a particular shape of C and C, shown on Fig. 1. The boundary conditions for the potential q5 are: I$ = t VI2
on A’B’
(6)
l$ = - v/2
on AB on AA’ and BB’.
(7)
vl#J*i&=o b(i) =
f
(8)
The analytical solution is obviously found to be:
It is also a classical problem of potential theory to 113
G. DE MEY
114
I
AC;
the algebraic set, a standard program was used based on the Gauss’ pivotal elimination method [20]. An estimation of the rank of the matrix could then also be performed. One will understand that if R becomes large, the offdiagonal elements of the matrix will attain the same order of magnitude as the diagonal elements. The matrix then becomes ill-conditioned. The influence on the solution is illustrated on Fig. 2, representing some numerical P
!,
a=b=lO.
I3
I3
0 R=?.25 q
+ R=lO. m R=12.5
1.
Fig. 1. Rectangular geometry with boundary C to solve the integral equation numerically. The curve C, has been chosen circular as the integral equation can then also be solved analytically.
o,~ 0.
determine the charge distribution p on the circle C, which causes a homogeneous field inside C,[18]. The
K/2
Fig. 2. Numerical and analytical results for the source function p (m = n = 10).
result is: results for p compared with the exact solution (10). If the dimensions of C, are about twice the dimensiohs of C, (IO) the error on p is already high. The fact that the errors for R = 7.25 are somewhat higher than the errors for R = 10 For the numerical solution of the problem, the curve is due to the fact that C, (R = 7.25) reaches almost the C, has been divided into n equal parts ACi (j = 1 . . n). corner points of C. If C, is close to C, the fundamental In each interval ACi the charge p is assumed to be a discontinuities for V~I should be taken into account, constant pi. Hence, the expression (4) for the potential 4 which has not been done here. Similar results are shown can be approximated by: in Table 1, representing Ap for various values of R where Ap is the largest difference between an exact and numerical p value normalised to the maximum p value. Table 1 gives also the “rank” of the matrix. If after the kth elimination step no pivot element can be found larger where !; is the centre point of AC; and [AC]the length of than E (= lo+‘) times the previous pivot element, the rank a segment. Imposing the boundary conditions (6)-(8) on of the matrix is said to be equal to k - 1. Because the pivot elements are never found to be exactly zero, the the proposed solution (ll), yields: numerical procedure can be continued until a result is obtained. However, the accuracy is low if the rank is less (174 than the dimension n. On Table 1, the error on the potential 4 is represented by means of the parameter A$J, which is the greatest difference between the exact and numerical +-value inside C after normalisation to the maximum value $J= 0.5. One concludes that the error Ad on the potential is much lower than the error Ap even for high values of the radius R. This is a very interesting feature of the method because the source function p is only an intermediate and not the final result of the problem. The low values In order to determine the II unknowns pi, we have to for 84 are explained by the fact that the potential 4 is place n points P, on the boundary C. These points are calculated by integrating p over the curve C,. The osalso placed at regular distances. The relations (12)-(14) cillatory behaviour of p (Fig. 2) will then be smoothed constitute then a linear set of algebraic equations which out, giving more accurate results of 4. For R > 10’ the can be easily solved numerically by the Gauss elimina- errors on p are no longer oscillatory so that large errors A4 are found (Table 1). tion method[19]. The potential 4 can then be evaluated Another advantage of using a curve CI different from by (11). The problem has been solved numerically for n = 40, C is the uniformity of the obtained accuracy. By using ordinary integral equation techniques (C = C,), the a = b and V = 1. The symmetry properties of the geometry have been taken into account in order to numerical discretisation process takes place on the bounreduce the matrix dimension to 10. For the solution of dary itself and the errors are larger in the neighborhood 2v p=-cos8. a
115
Integral equations for potential problems Table 1. R
AP
rank
1.25 10 12.5 15 17.5 20 25 30 40 50 75 100 160 200 500 103 3.103 lo4 3.104 105 3.105 lo6 3.106
0.1 10-z 0.5 2 2.5 0.5 10 43 2.5 5.1 7.8 6.2 12 1.3 19 0.8 1.5 0.8 5.2 4.7 14 0.5 1.6
10 10 10 9 7 8 7 8 7 6 7 6 6 6 6 6 6 5 5 5 5 5 4
A4 0.2 2.10-5 10-h 10-e 2.10-6 1.5 10-h 5.10-6 5.10-5 2.10-e 10-x 10-5 10-4 1.2 10-d 5.10-5 5.10-4 7.10-S 6.10-3 3.10-3 6.10-* 7.10-2 0.6 0.375 1 ?< I.LJ
of the boundary. By the method outlined here, the accuracy does not vary appreciably inside C. Even for points on C, the potential do can be calculated with the same accuracy. 3. CONCLUSION In this paper a numerical method has been presented
to solve Laplace equation in a two dimensional area. The method is based on an integral equation technique where the source function is not defined on the boundary C but on another curve C, surrounding C. In spite of the fact that the matrix of the algebraic set can be ill conditioned, the accuracy for the calculated potential is very good, even when the numerically calculated source function shows large errors. Due to the ill-conditioned matrix a pivotal elimination method is necessary. This procedure allows to get an estimation of the rank of the matrix. From the example treated above, one can conclude that the dimensions of C, may be up to 1000times larger than the dimensions of C in order to get a sufficient accuracy. The accuracy improves if C, is closer to C. However, if CI reaches C too closely, the error still increases because the integral equation is not corrected for the fundamental discontinuities. By comparing this method
with the usual integral equation method (C = C,), it is found that the accuracy is uniformly distributed over the area inside C. REFERENCES
1. J. Van Blade], Electromagnetic Fields pp. 131-133McGraw Hill, New York (1964). 2. 3. Van Blade], Low frequency scattering through an aperture in a rigid screen. J. Sound Vibration 6, 386-395 (1967). 3. J. Van Blade], Low frequency scattering through an aperture in a soft screen. J. Sound Vibration 8, 186-195 (1968). 4. H. J. Pauwels and G. De Mey, Surface to surface transition probabilities in thin filmcapacitors. Physica Status Solidi. a 24, K39-K44 (1974). 5. 0. Zienkiewicz, The Finite Element Method in Engineering Science p. 272 McGraw Hill, London (1971). 6. T. W. Edwards and J. Van Bladel, Electrostatic dipole moment of a dielectric cube. Appl. Sci. Res. 9, 151-155 (1961). 7. K. Mei and J. Van Blade], Low frequency scattering by rectangular cylinders. IEEE Truns. Antennas Propagation AP-11, 52-56 (1963). 8. K. Mei and J. Van Blade], Scattering by perfectly conducting rectangular cylinders. IEEE Trans. Antennas Propagation AP-11, 185-192(1963). 9. R. Shaw. An integral equation approach to diffusion. Inl. J. Heat Mass Transfer 17, 693-699 (1974).
10. G. De Mey, A method for calculating eddy currents in plates of arbitrary geometry. Arch. Elektrotkchn. 56, 137-140 (1974). 11. G. De Mey, Frequency shift in resonant cavities. Int. J. Electron. 37, 369-315 (1974).
12. G. De Mey, Integral equation for the potential distribution in a Hall generator. Ekctronics L&t. 9, 264-266 (1973). 13. G. De Mey, Determination of the electric field in a Hall generator under influence of an alternating magnetic field. Solid State Electron. 17, 977-919 (1974).
14. G. De Mey, Calculation of eigenvalues of the Helmholtz equation by an integral equation. Numer. Merh. Engng 10, 59-66 (1976). 15. G. De Mey, An integral equation method to solve Poisson’s equation for stochastic charge densities and boundary potential. Lett. Nuovo Cimento 12, 47U72 (1975). 16. G. De Mey, An integral equation method for the numerical calculation of ion drift and diffusion in evaporated dielectrics. Computing 17, 169-176 (1976). 17. J. Van Blade], Electromagnetic Fields p. 57-59 McGraw Hill, New York (1%4). 18. W. T. Scott, The Physics of Electricity and Magnetism. pp. 168-170Wiley, New York (1962). 19. T. R. McCalIa, Introduction to Numerical methods and FORTRAN Programming p. 171 Wiley, New York (1967). 20. I.B.M.: System 360 Scientific subroutine package. p, 121-124 IBM Technical Publications, New York, (1970). 21. S. De Wolf and G. De Mey, Numerical methods for solving integral equations of potential problems. Inform. Proc. Lett. 3, 121-124(1975).
Vol. 8. pp. ll?-115.
Pergamon Press 1978.
Printed in Great Britain
INTEGRAL EQUATIONS FOR POTENTIAL PROBLEMS WITH THE SOURCE FUNCTION NOT LOCATED ON THE BOUNDARY G. DE MEY
Laboratory of Electronics, Ghent State University, Sint Pietersnieuwstraat 41, 9000 Ghent, Belgium (Received December 1975)
Abstract-A special integral equation method for solving Laplace equation in a two dimensional area is studied. Usually an unknown source function is defined along the boundary. In the method to be outlined here, this source function is located on another closed curve. This new method offers several advantages as compared to the classical integral equation technique.
1.INTRODUCTION
Many problems in physics and engineering which are usually described by partial differential equations and appropriate boundary conditions can also be described by equivalent integral equations [l-5]. Integral equations have proved to be very useful for the numerical solution of a given problem&161. If the original problem is two dimensional, the corresponding integral equation will be one-dimensional. This fact reduces the complexity of the computer programming and the computation time is also diminished. For the particular case of a two dimensional potential problem described by the Laplace equation: V’4 = 0
(1)
in an area S with boundary C, the potential 4 is written as: 4(f) = f p(i’)G($‘) dC
(2)
where G(ili’)=llnL
2Tr JJ- F”)
is the Green’s function of the eqn (l)[l]. p(fi is an unknown source function defined along the boundary C. By imposing the boundary conditions on the proposed solution (3), one obtains an integral equation for the function p(fl. Note that the Laplace eqn (1) has been replaced by a one-dimensional integral equation. In the expression (2), the function p(F) is defined along the boundary C of S. However, there is no a prior’ reason why the function p(!) should be defined along C. Taking a second closed line C,, which completely surrounds the area S, one can write the potential 4 as:
the unknown function p(Q, which is now defined along C,. If C, does not coincide with C, the fundamental discontinuities for a boundary condition such as V&E” = 0 should not be taken into account[l7]. One can simply write that: V4(n =
f
p(i')G(i/i') dC;.
(4)
Because the Green’s function (3) is a particular solution of (I), the expression (4) is also a solution of Laplace equation. We have now to impose the boundary condition on (4), in order to obtain an integral equation for CAS Vol. 8, No. 1-H
p(Y)VG(@‘) dC;
(5)
where (5) is still valid on the boundary C. Another feature caused by introducing C, is related to the numerical solution of the integral equation. The continuous function p(?) is then replaced by a set of n numbers pI . . . p. along C,. If C, coincides with C, the numerical noise will be large in the vicinity of the boundary, because the numerical discretisation process takes place there. If C, does not coincide with C, the numerical noise can be made sufficiently low provided C, is not too close to C. This problem will be investigated in this paper for a particular example where C1 is a circular and C a rectangular boundary shape. The boundary conditions on C are chosen in such a way that the problem can also be solved analytically. The accuracy of the numerical results can then easily be checked. 2. ExAMF%E
The numerical solution of Laplace equation in an area S with boundary C by means of an integral equation with the source function defined along a closed line C, (Z C) will be given for a particular shape of C and C, shown on Fig. 1. The boundary conditions for the potential q5 are: I$ = t VI2
on A’B’
(6)
l$ = - v/2
on AB on AA’ and BB’.
(7)
vl#J*i&=o b(i) =
f
(8)
The analytical solution is obviously found to be:
It is also a classical problem of potential theory to 113
G. DE MEY
114
I
AC;
the algebraic set, a standard program was used based on the Gauss’ pivotal elimination method [20]. An estimation of the rank of the matrix could then also be performed. One will understand that if R becomes large, the offdiagonal elements of the matrix will attain the same order of magnitude as the diagonal elements. The matrix then becomes ill-conditioned. The influence on the solution is illustrated on Fig. 2, representing some numerical P
!,
a=b=lO.
I3
I3
0 R=?.25 q
+ R=lO. m R=12.5
1.
Fig. 1. Rectangular geometry with boundary C to solve the integral equation numerically. The curve C, has been chosen circular as the integral equation can then also be solved analytically.
o,~ 0.
determine the charge distribution p on the circle C, which causes a homogeneous field inside C,[18]. The
K/2
Fig. 2. Numerical and analytical results for the source function p (m = n = 10).
result is: results for p compared with the exact solution (10). If the dimensions of C, are about twice the dimensiohs of C, (IO) the error on p is already high. The fact that the errors for R = 7.25 are somewhat higher than the errors for R = 10 For the numerical solution of the problem, the curve is due to the fact that C, (R = 7.25) reaches almost the C, has been divided into n equal parts ACi (j = 1 . . n). corner points of C. If C, is close to C, the fundamental In each interval ACi the charge p is assumed to be a discontinuities for V~I should be taken into account, constant pi. Hence, the expression (4) for the potential 4 which has not been done here. Similar results are shown can be approximated by: in Table 1, representing Ap for various values of R where Ap is the largest difference between an exact and numerical p value normalised to the maximum p value. Table 1 gives also the “rank” of the matrix. If after the kth elimination step no pivot element can be found larger where !; is the centre point of AC; and [AC]the length of than E (= lo+‘) times the previous pivot element, the rank a segment. Imposing the boundary conditions (6)-(8) on of the matrix is said to be equal to k - 1. Because the pivot elements are never found to be exactly zero, the the proposed solution (ll), yields: numerical procedure can be continued until a result is obtained. However, the accuracy is low if the rank is less (174 than the dimension n. On Table 1, the error on the potential 4 is represented by means of the parameter A$J, which is the greatest difference between the exact and numerical +-value inside C after normalisation to the maximum value $J= 0.5. One concludes that the error Ad on the potential is much lower than the error Ap even for high values of the radius R. This is a very interesting feature of the method because the source function p is only an intermediate and not the final result of the problem. The low values In order to determine the II unknowns pi, we have to for 84 are explained by the fact that the potential 4 is place n points P, on the boundary C. These points are calculated by integrating p over the curve C,. The osalso placed at regular distances. The relations (12)-(14) cillatory behaviour of p (Fig. 2) will then be smoothed constitute then a linear set of algebraic equations which out, giving more accurate results of 4. For R > 10’ the can be easily solved numerically by the Gauss elimina- errors on p are no longer oscillatory so that large errors A4 are found (Table 1). tion method[19]. The potential 4 can then be evaluated Another advantage of using a curve CI different from by (11). The problem has been solved numerically for n = 40, C is the uniformity of the obtained accuracy. By using ordinary integral equation techniques (C = C,), the a = b and V = 1. The symmetry properties of the geometry have been taken into account in order to numerical discretisation process takes place on the bounreduce the matrix dimension to 10. For the solution of dary itself and the errors are larger in the neighborhood 2v p=-cos8. a
115
Integral equations for potential problems Table 1. R
AP
rank
1.25 10 12.5 15 17.5 20 25 30 40 50 75 100 160 200 500 103 3.103 lo4 3.104 105 3.105 lo6 3.106
0.1 10-z 0.5 2 2.5 0.5 10 43 2.5 5.1 7.8 6.2 12 1.3 19 0.8 1.5 0.8 5.2 4.7 14 0.5 1.6
10 10 10 9 7 8 7 8 7 6 7 6 6 6 6 6 6 5 5 5 5 5 4
A4 0.2 2.10-5 10-h 10-e 2.10-6 1.5 10-h 5.10-6 5.10-5 2.10-e 10-x 10-5 10-4 1.2 10-d 5.10-5 5.10-4 7.10-S 6.10-3 3.10-3 6.10-* 7.10-2 0.6 0.375 1 ?< I.LJ
of the boundary. By the method outlined here, the accuracy does not vary appreciably inside C. Even for points on C, the potential do can be calculated with the same accuracy. 3. CONCLUSION In this paper a numerical method has been presented
to solve Laplace equation in a two dimensional area. The method is based on an integral equation technique where the source function is not defined on the boundary C but on another curve C, surrounding C. In spite of the fact that the matrix of the algebraic set can be ill conditioned, the accuracy for the calculated potential is very good, even when the numerically calculated source function shows large errors. Due to the ill-conditioned matrix a pivotal elimination method is necessary. This procedure allows to get an estimation of the rank of the matrix. From the example treated above, one can conclude that the dimensions of C, may be up to 1000times larger than the dimensions of C in order to get a sufficient accuracy. The accuracy improves if C, is closer to C. However, if CI reaches C too closely, the error still increases because the integral equation is not corrected for the fundamental discontinuities. By comparing this method
with the usual integral equation method (C = C,), it is found that the accuracy is uniformly distributed over the area inside C. REFERENCES
1. J. Van Blade], Electromagnetic Fields pp. 131-133McGraw Hill, New York (1964). 2. 3. Van Blade], Low frequency scattering through an aperture in a rigid screen. J. Sound Vibration 6, 386-395 (1967). 3. J. Van Blade], Low frequency scattering through an aperture in a soft screen. J. Sound Vibration 8, 186-195 (1968). 4. H. J. Pauwels and G. De Mey, Surface to surface transition probabilities in thin filmcapacitors. Physica Status Solidi. a 24, K39-K44 (1974). 5. 0. Zienkiewicz, The Finite Element Method in Engineering Science p. 272 McGraw Hill, London (1971). 6. T. W. Edwards and J. Van Bladel, Electrostatic dipole moment of a dielectric cube. Appl. Sci. Res. 9, 151-155 (1961). 7. K. Mei and J. Van Blade], Low frequency scattering by rectangular cylinders. IEEE Truns. Antennas Propagation AP-11, 52-56 (1963). 8. K. Mei and J. Van Blade], Scattering by perfectly conducting rectangular cylinders. IEEE Trans. Antennas Propagation AP-11, 185-192(1963). 9. R. Shaw. An integral equation approach to diffusion. Inl. J. Heat Mass Transfer 17, 693-699 (1974).
10. G. De Mey, A method for calculating eddy currents in plates of arbitrary geometry. Arch. Elektrotkchn. 56, 137-140 (1974). 11. G. De Mey, Frequency shift in resonant cavities. Int. J. Electron. 37, 369-315 (1974).
12. G. De Mey, Integral equation for the potential distribution in a Hall generator. Ekctronics L&t. 9, 264-266 (1973). 13. G. De Mey, Determination of the electric field in a Hall generator under influence of an alternating magnetic field. Solid State Electron. 17, 977-919 (1974).
14. G. De Mey, Calculation of eigenvalues of the Helmholtz equation by an integral equation. Numer. Merh. Engng 10, 59-66 (1976). 15. G. De Mey, An integral equation method to solve Poisson’s equation for stochastic charge densities and boundary potential. Lett. Nuovo Cimento 12, 47U72 (1975). 16. G. De Mey, An integral equation method for the numerical calculation of ion drift and diffusion in evaporated dielectrics. Computing 17, 169-176 (1976). 17. J. Van Blade], Electromagnetic Fields p. 57-59 McGraw Hill, New York (1%4). 18. W. T. Scott, The Physics of Electricity and Magnetism. pp. 168-170Wiley, New York (1962). 19. T. R. McCalIa, Introduction to Numerical methods and FORTRAN Programming p. 171 Wiley, New York (1967). 20. I.B.M.: System 360 Scientific subroutine package. p, 121-124 IBM Technical Publications, New York, (1970). 21. S. De Wolf and G. De Mey, Numerical methods for solving integral equations of potential problems. Inform. Proc. Lett. 3, 121-124(1975).