- Email: [email protected]
Numerical methods for solving integral equations of potential problems
Volume 3, number 4
INFORMATION BPOCEWNG
LETTERS
March,1975
a
ERIC
HODS FOR SOLVING INTEGRAL
WOJ,F and G. De MEY rut09 of Ekectrvnics, Ghent State Urdversity, St-Rctersnieuwstraat 41, B-9000 Ghent, Belgium Received 2 December 1974
potential problems
integral equations
roblems in physics and engineering which described by partial differential equations can also be treated by equivalent integral equations. We cite e.g. electrostatic potential [ 1,2], acoustic [3,4] s electromagnetic scattering [S, 61, elasticity diffusion phenomena [8 ,9 ] and eigeqvalue problems [ 10,l J]. By using an integral equation, the dimension of the problem is reduced by one. This will save computati~ timeif the problem is to be calculated numeriCaJlYIn this paper an integral equation for s specifk potential problem will be solved for different shapes of the unknown boundary function. ‘Jhe \ . .;ious results will be compared on the.basis of their accuracy.
ii
nt
We consider a two-dimensional electrostatic potential problem in the area S with boundary C (fig. 1):
v24$=0.
(0
The metallic boundaries Z1 and Cj are held at constant component of 0 e no 5 should vanish at the other boundaries C2 and Cd :
Fig. 1. Rectangular geometry used to compare the numerical results with the analytical solution.
Using the Green’s funcilcm
lNFORMA’ITON PROCESSING LETTJZKS
Vdw#u? 3, number 4
)
=fp(s’)G(slr’W’ ,
(6)
C
re p(r) is a function along the boundary C. Impothe boundary conditions (2), (S) and (4) on the d solution (6), gives result to: 1 vJ=&,
PV)
In ._-l-_ jr _ #l CiC’,
roNI;
(7)
f
‘ionC? andC4.
W
(71, (8) and (9 5ccnstitute an integral equation for ;le unknown function p(r). The first term appearilt (9) is caused by the &continuity of VqP tin [ 121. p(t) kr been determine& the potential @(Y)can culated by using formula (6). e te&iques which will be outlined in this paper can be ap$ied to arbitrary shapes of the boundary C. flowever, we shall restrict ourselves to a simpte rectr geometpr (fig. 1) in order to compare numericts with thf:;analytical solution.
Fig. 2. Difference A between the numericaUy and anratytiay calculated potentiai for VI = -I5 = IO as obtained by the four methods & B, C and D. Each side of the square geomb try has been divided into four segments. Curve a represents A in the points (0.9a + 0.1 Ro, -0.9b) for A= 1, .... 9. Curve b represents A in the points (0.94, -0.9b + 0.1 kb) for k = 1, -**9, l
In order to solve the integral equation (7)-(g) numerically, the boundary C is divided into n segments ACj. The continuous function p(i) is then replaced n constants PI, ... . pn. The integral equation becorn&san n I( n linear algebraic set. Four methods, dencrted furtheron as A, B, C:and D will be discussed arm compared on the basis of their accuracy. 1 The flJst method A has already been used to determine the potential distribution in a Hall generator (13). The methods A and B are now explained in detail. In each segment ACj the firnction p is assumed 10be a confitant pi, hence the exnression (6) for the Frtential # can be written as:
The integral equation can be treated in a similar fashion. For (7), we get: VI
=kjq
Pj
O”
SC(rilr’)dC’,
ri onC1.
(11)
AC* I
(11) should be fuifiUed for all points Q which are the centre points of the segments lying on the boundary Cl. (8) and (9) are treated similarly. We have now obtained a linear algebraic set for the pi unknowns ~1, .... pn. For the method B, the integrals appearing in (10) and (11) are calculated exactly, the remaining terms being approximated by
where jACjI denotes &e length of the jth segment ‘142
INFORMATION PROCESSING LETTERS
March 19?5
pointsYi are re9.L 6 for the ;lumerical solution 4
(AJ?tiW?ARV WITS>
o:F
the integral e+
:I- i i I)* Time points have been chosen & the end points of the e$ments as l‘arthey do not coincide with one of the comers. For these segments,new check points Year placed at 0.3 IAC’I from the corr9rs. The resuhs are more accurate than those obtained with the other methods (fig. 2) and the oscillations have disappeared (fig. 3).
4. Conclusion
FJig.3. D for a
density p calculated by the methoda A, B, C and metry with a/b = S/7. Note the oscihtionrin the c& f&apondiag to C.
A$ The diagonal terms cannot be calculatedby (12), due to the singularityof the Green’sfunction (5) if q=q.For Vt = - Vz = 10, the problemhas been solved for II= 16 (4 segments on each side) by the Gauss pivotal e?Wnation method [ 141. Figs. 2a and 2b show the difference A between the exact solution lOy/b and the numerical resultsalong the lines x = 0.9a en y = 0.9b respectively,as the errorsare largestclose to the boundary. Increasingn leads also to a remarkablerise of accuracy. Results for the function p are shown in fig. 3. For the method C, p is assumed to vzuylinearly from p to pel in each segmenl:ACj (f%. 1). At it is believed &at this continuous and piecewise linear function would be a better approximation,more accurateresults afe expected. However,numericalresults for p show an oscillatory behaviour(fig. 3), which is certainly incorrect. Neverthelessthe erroron the potential @is still of the same order of magnitude as the resultsobtained with tis:emethods A and B (fig. 2). It has been found that the oscillations are caused by the built-in continuity of p even at a cornerpoint, where the boundary condition changes. In method D, the function p is approximatedby a piecewise linear function in each segment but a discontinuity at a comer point is no 1ongc:rexcluded. At a cornerpoint two values pj and ,ei sh&d then be determined.For the situation depicted uafig. 1, xhe total number of unknowns is then n + 4. for n segments. Even rr + 4 check
As a general conclusion we may state that the accuracy increases if p is approximated by a piecewise linear fu action (methcd D) in stead of by constant values (method R) or Dirac pulses for the non-diagonal tern-s (method A). Similar results have been obtained with the fmite element method where tie . accuracy increases if the trial function is 8 polynomial of higher degree [ 151. However, one Grould take into account that the unknown function p can undergo a discontinuity at a point where the boundary condition changes. Otherwise, oscillatory solutions will be obtained (method C), Nevertheless, the accuracy for the potential 0 obtained with method C is quite good, which proves that possible oscillations in the function p will be filtered out due to the integration (10). This particular effect makes the integral equation more attractive than other numerical methods based on finite difference approximations where possible instabilities or numerical noise will not be smoothed by integrations [ 161.
Acknowledgements The authors wish to thank the IWONL of Belgium, who has partially sponsored this project.
References (11T.W. Edwards and J. Van Blade&Electrostatic dipole moment of a dielectric cube, Appl. SC. Res. 9 (1961) J5 l155. VI G. De Mey, A method for calculating eddy currents in plates of arbitrary geometry, Archiv far Elektrotechnik 56 (1974) 137-140. 831 J. Van Blade&Low frequency scattering through an aperr. 6 (1967) 386-395. screen, J. Soun ture in a
123
V&n&e 3, number 4
INFORMATION PROCESSINGLETTERS
1[41 1. Van Blade&Low frequency scattering through an aperture in a soft screen, J. Sound and Vibr. 8 t 1968) 186195. 151 K. Md and Z. Van Blade&Low frequency scattering by cylinders, lEEE Trans. Antennar’Prop., 3) 52.-56. 16) K. Mei and J. Van Blade&Scattering by perfectly condu&ing rectangular cylinders, IEEE Trans. Antennas Rroy. AP-I l(1963) 185-192. , The finite element method in engineering w Iii& London, 1971) 272. egraf equation approach to diffusion, Intern. 1. WeatMassTransfer 17 (1974) 693-699. 19) N.J. R uwels and G. De Mey, Surface to surface transitIotrprobabilities in thin film capacitr..rs,JWjs. Stat. Sol. 3 24 (1974) K J9--K44, 110) G. De Mey, Frequency shift in r!:sonant cavities, Interu. J. of EJcctronics 37 (1974) 369--375.
March 1975
WI G. De Mey, Calculation of eigenvaluts of the Helmholtz equation by an integral, Numerical methods in Engiaeering (to be published). 1121J. Van Bladel, Electromagnetic fields, MCCraw HU, New York (1964) p. 58. 1131 C. De Mey, Integral equation for the potential distribution in a Hall generator, Electronics Letters 9 (1973) 264-266. 114) I.B.M, Application Program: System 360 rdentlfic subroutine package, version 3 (I.B.M. Technical Publications, New York, 1970) p. 121-124. 0. Zienkiewicx, The finite element method in engineering science (McGraw HipI,London,. 1971) 103-120. WI G.D. Smith: Numerical solution of partial differential equations (Oxford University Press, London, 1971) 132-140.
INFORMATION BPOCEWNG
LETTERS
March,1975
a
ERIC
HODS FOR SOLVING INTEGRAL
WOJ,F and G. De MEY rut09 of Ekectrvnics, Ghent State Urdversity, St-Rctersnieuwstraat 41, B-9000 Ghent, Belgium Received 2 December 1974
potential problems
integral equations
roblems in physics and engineering which described by partial differential equations can also be treated by equivalent integral equations. We cite e.g. electrostatic potential [ 1,2], acoustic [3,4] s electromagnetic scattering [S, 61, elasticity diffusion phenomena [8 ,9 ] and eigeqvalue problems [ 10,l J]. By using an integral equation, the dimension of the problem is reduced by one. This will save computati~ timeif the problem is to be calculated numeriCaJlYIn this paper an integral equation for s specifk potential problem will be solved for different shapes of the unknown boundary function. ‘Jhe \ . .;ious results will be compared on the.basis of their accuracy.
ii
nt
We consider a two-dimensional electrostatic potential problem in the area S with boundary C (fig. 1):
v24$=0.
(0
The metallic boundaries Z1 and Cj are held at constant component of 0 e no 5 should vanish at the other boundaries C2 and Cd :
Fig. 1. Rectangular geometry used to compare the numerical results with the analytical solution.
Using the Green’s funcilcm
lNFORMA’ITON PROCESSING LETTJZKS
Vdw#u? 3, number 4
)
=fp(s’)G(slr’W’ ,
(6)
C
re p(r) is a function along the boundary C. Impothe boundary conditions (2), (S) and (4) on the d solution (6), gives result to: 1 vJ=&,
PV)
In ._-l-_ jr _ #l CiC’,
roNI;
(7)
f
‘ionC? andC4.
W
(71, (8) and (9 5ccnstitute an integral equation for ;le unknown function p(r). The first term appearilt (9) is caused by the &continuity of VqP tin [ 121. p(t) kr been determine& the potential @(Y)can culated by using formula (6). e te&iques which will be outlined in this paper can be ap$ied to arbitrary shapes of the boundary C. flowever, we shall restrict ourselves to a simpte rectr geometpr (fig. 1) in order to compare numericts with thf:;analytical solution.
Fig. 2. Difference A between the numericaUy and anratytiay calculated potentiai for VI = -I5 = IO as obtained by the four methods & B, C and D. Each side of the square geomb try has been divided into four segments. Curve a represents A in the points (0.9a + 0.1 Ro, -0.9b) for A= 1, .... 9. Curve b represents A in the points (0.94, -0.9b + 0.1 kb) for k = 1, -**9, l
In order to solve the integral equation (7)-(g) numerically, the boundary C is divided into n segments ACj. The continuous function p(i) is then replaced n constants PI, ... . pn. The integral equation becorn&san n I( n linear algebraic set. Four methods, dencrted furtheron as A, B, C:and D will be discussed arm compared on the basis of their accuracy. 1 The flJst method A has already been used to determine the potential distribution in a Hall generator (13). The methods A and B are now explained in detail. In each segment ACj the firnction p is assumed 10be a confitant pi, hence the exnression (6) for the Frtential # can be written as:
The integral equation can be treated in a similar fashion. For (7), we get: VI
=kjq
Pj
O”
SC(rilr’)dC’,
ri onC1.
(11)
AC* I
(11) should be fuifiUed for all points Q which are the centre points of the segments lying on the boundary Cl. (8) and (9) are treated similarly. We have now obtained a linear algebraic set for the pi unknowns ~1, .... pn. For the method B, the integrals appearing in (10) and (11) are calculated exactly, the remaining terms being approximated by
where jACjI denotes &e length of the jth segment ‘142
INFORMATION PROCESSING LETTERS
March 19?5
pointsYi are re9.L 6 for the ;lumerical solution 4
(AJ?tiW?ARV WITS>
o:F
the integral e+
:I- i i I)* Time points have been chosen & the end points of the e$ments as l‘arthey do not coincide with one of the comers. For these segments,new check points Year placed at 0.3 IAC’I from the corr9rs. The resuhs are more accurate than those obtained with the other methods (fig. 2) and the oscillations have disappeared (fig. 3).
4. Conclusion
FJig.3. D for a
density p calculated by the methoda A, B, C and metry with a/b = S/7. Note the oscihtionrin the c& f&apondiag to C.
A$ The diagonal terms cannot be calculatedby (12), due to the singularityof the Green’sfunction (5) if q=q.For Vt = - Vz = 10, the problemhas been solved for II= 16 (4 segments on each side) by the Gauss pivotal e?Wnation method [ 141. Figs. 2a and 2b show the difference A between the exact solution lOy/b and the numerical resultsalong the lines x = 0.9a en y = 0.9b respectively,as the errorsare largestclose to the boundary. Increasingn leads also to a remarkablerise of accuracy. Results for the function p are shown in fig. 3. For the method C, p is assumed to vzuylinearly from p to pel in each segmenl:ACj (f%. 1). At it is believed &at this continuous and piecewise linear function would be a better approximation,more accurateresults afe expected. However,numericalresults for p show an oscillatory behaviour(fig. 3), which is certainly incorrect. Neverthelessthe erroron the potential @is still of the same order of magnitude as the resultsobtained with tis:emethods A and B (fig. 2). It has been found that the oscillations are caused by the built-in continuity of p even at a cornerpoint, where the boundary condition changes. In method D, the function p is approximatedby a piecewise linear function in each segment but a discontinuity at a comer point is no 1ongc:rexcluded. At a cornerpoint two values pj and ,ei sh&d then be determined.For the situation depicted uafig. 1, xhe total number of unknowns is then n + 4. for n segments. Even rr + 4 check
As a general conclusion we may state that the accuracy increases if p is approximated by a piecewise linear fu action (methcd D) in stead of by constant values (method R) or Dirac pulses for the non-diagonal tern-s (method A). Similar results have been obtained with the fmite element method where tie . accuracy increases if the trial function is 8 polynomial of higher degree [ 151. However, one Grould take into account that the unknown function p can undergo a discontinuity at a point where the boundary condition changes. Otherwise, oscillatory solutions will be obtained (method C), Nevertheless, the accuracy for the potential 0 obtained with method C is quite good, which proves that possible oscillations in the function p will be filtered out due to the integration (10). This particular effect makes the integral equation more attractive than other numerical methods based on finite difference approximations where possible instabilities or numerical noise will not be smoothed by integrations [ 161.
Acknowledgements The authors wish to thank the IWONL of Belgium, who has partially sponsored this project.
References (11T.W. Edwards and J. Van Blade&Electrostatic dipole moment of a dielectric cube, Appl. SC. Res. 9 (1961) J5 l155. VI G. De Mey, A method for calculating eddy currents in plates of arbitrary geometry, Archiv far Elektrotechnik 56 (1974) 137-140. 831 J. Van Blade&Low frequency scattering through an aperr. 6 (1967) 386-395. screen, J. Soun ture in a
123
V&n&e 3, number 4
INFORMATION PROCESSINGLETTERS
1[41 1. Van Blade&Low frequency scattering through an aperture in a soft screen, J. Sound and Vibr. 8 t 1968) 186195. 151 K. Md and Z. Van Blade&Low frequency scattering by cylinders, lEEE Trans. Antennar’Prop., 3) 52.-56. 16) K. Mei and J. Van Blade&Scattering by perfectly condu&ing rectangular cylinders, IEEE Trans. Antennas Rroy. AP-I l(1963) 185-192. , The finite element method in engineering w Iii& London, 1971) 272. egraf equation approach to diffusion, Intern. 1. WeatMassTransfer 17 (1974) 693-699. 19) N.J. R uwels and G. De Mey, Surface to surface transitIotrprobabilities in thin film capacitr..rs,JWjs. Stat. Sol. 3 24 (1974) K J9--K44, 110) G. De Mey, Frequency shift in r!:sonant cavities, Interu. J. of EJcctronics 37 (1974) 369--375.
March 1975
WI G. De Mey, Calculation of eigenvaluts of the Helmholtz equation by an integral, Numerical methods in Engiaeering (to be published). 1121J. Van Bladel, Electromagnetic fields, MCCraw HU, New York (1964) p. 58. 1131 C. De Mey, Integral equation for the potential distribution in a Hall generator, Electronics Letters 9 (1973) 264-266. 114) I.B.M, Application Program: System 360 rdentlfic subroutine package, version 3 (I.B.M. Technical Publications, New York, 1970) p. 121-124. 0. Zienkiewicx, The finite element method in engineering science (McGraw HipI,London,. 1971) 103-120. WI G.D. Smith: Numerical solution of partial differential equations (Oxford University Press, London, 1971) 132-140.