Matlab about Numerical Solution of Laplace Equation

1 Numerical Solution of Laplace Equation The electric potential in a charge-free region is given by the Laplace equation ∇2 V = 0 (1) In Cartesian coordinates in two dimensions, this equation expands a s 2 V x 2+ 2 V y 2= 0 (2) An analytical solution of this partial di erential equation can be fou nd only if the domain of interest (i.e., the region in which the solution is to be computed) is simple, such as a rectangle or a circle. For nding the vo ltage distribution in a region of nonuniform geometry, engineers frequen tly use numerical techniques.

There are several numerical methods that can be used for solving equation (1) among which the nite di erence method method is the simplest. T he following is a brief introduction to the nite di erence method.

The nite di erence method consists of three steps:

• Discretize the domain as a grid.

• Find an algebraic equation for each node in the grid.

• Solve the set of algebraic equations. 1 Grid2 1 Grid For nite di erence solution, we direct our attention only at some nite number of points in the domain. These points are referred to as grid (or mesh or node) points. For domains in Cartesian geometry, a rectan gular grid is usually used. For cylindrical geometry, the region is discretized ut ilizing cylindrical coordinates. x y Fig 1. Discretization It is not necessary that the grid size be uniform for the xand y coordinates, however a uniform grid results in some simpli cation in th e numerical method. Note however that for domains with irregular sh aped geometry, a uniform grid cannot be used. Nevertheless, it is possib le to extend the method with good accuracy.

The potential Vat the boundary nodes are given. The objective is to nd the potential Vfor the interior nodes. 2 Di erence Equation3 2 Di erence Equation This step generates an algebraic equation for each node in the domain. This is based on the Taylor series expansion:

V(x + x) = V(x ) + V x( x ) x+ 1 2 2 V x 2( x )( x)2 + higher order terms V (x − x) = V(x ) − V x ( x ) x+ 1 2 2 V x 2( x )( x)2 + higher order terms (3) Adding the two equations and neglecting higher order terms, we obt ain 2 V x 2( x ) = V (x + x) − 2V (x ) + V(x − x) ( x)2 (4) Note that Vis also a function of yso that the above equation actually holds for any given y. Thus we rewrite the above equation as 2 V x 2( x, y ) = V (x + x, y)− 2V (x, y ) + V(x − x, y ) ( x)2 (5) Similar expansion can be done for the potential in the y-coordinate to obtain 2 V y 2( x, y ) = V (x, y + y) − 2V (x, y ) + V(x, y − y) ( y)2 (6) Substituting the above two equations in (2), we obtain V (x + x, y)− 2V (x, y ) + V(x − x, y ) ( x)2 +V (x, y + y) − 2V (x, y ) + V(x, y − y) ( y)2 = 0 (7) If we choose x= y, the above equation simpli es to V (x + x, y) +V(x − x, y ) +V(x, y + y) + V(x, y − y)− 4V (x, y ) = 0 (8) Denote the grid points in xcoordinate as i= 1 ,2 , , N xand those for the y-coordinate as j= 1 ,2 , , N y. Then the above equation is expressed as Vi+1 ,j+ V i 1,j + V i,j +1 + V i,j 1 − 4V i,j = 0 (9) 2 Di erence Equation4 x y i-1,j i+1,j i,j+1 i,j-1 i,j Fig 2. Di erence Formula This equation is called the di erence equation that approximates the Laplace equation at the grid point ( i, j). Clearly it shows that the potential V at any grid point is the average of that of four neighboring nodes.

For the numerical solution of the Laplace equation, a di erence equ ation must be derived for each interior node of the domain which are then s olved simultaneously for the unknown potential V.

Note also that in deriving the di erence formula (9), we have neglect ed all higher order terms in the Taylor series expansion. As such, the numerical solution obtained by this method is only an approximate solu tion.

Nevertheless, accuracy of the computed solution can be improved by discretizing the domain to a ner grid, i.e. by taking xand ysmall.

This also means that there will be many more equations that must be s olved simultaneously which increases computer time and demands more mem ory. 3 Solution5 3 Solution The next step is to solve the di erence equations derived in the previous step.

Again, there are various ways these equations can be solved. For e xample, one may write the di erence equations as a linear matrix equation and nd the solution by inverting the matrix equation. An easier alternative is to nd the solution using an iterative method.

For the iterative solution, rst assume some values for the unknow n potential Vfor the interior nodes; in fact V i,j = 0 is also a possible initial guess of the solution. Then sequentially update the solution for eac h node by taking the average of four nearest neighboring nodes:

Vi,j = 1 4 V i+1 ,j+ V i 1,j + V i,j +1 + V i,j 1 (10) Assume a rectangular domain with potential on its four boundaries a re known. This means that potential at all boundary nodes are known , and we have to calculate the potential at interior nodes only. Based on the above analysis, here is the pseudo-code for the algorithm:

1. Assume initial guess of the solution 2. Update the solution using for i=2 to N x− 1 for j=2 to N y− 1 V (i, j ) = 1 4 V i+1 ,j+ V i 1,j + V i,j +1 + V i,j 1 end for end for 3. Repeat from Step 2 until the solution converges For convergence test, take the sum of square of potential of all nodes, i.e., i jV 2 i,j for each iteration, and compare with that computed in the previous iteration. Stop the iteration when increment of the squar ed sum of potentials for two successive iterations is small (say 10 5 ).

To calculate the electric eld intensity, we use the equation E= ∇ V (11) 3 Solution6 For the two dimensional geometry, this equation simpli es to E= − V x a x − V y a y (12) Using equations (3), we then obtain E (x, y ) = −V (x + x, y)− V(x − x, y ) 2 x a x− V (x, y + y) − V(x, y − y) 2 y a y (13) A pseudo-code for computer implementation of this calculation is as f ollows for i=2 to N x− 1 for j=2 to N y− 1 E (i, j ) = [ −V i+1 ,j V i 1,j 2 x , − V i,j +1 V i,j 1 2 y ] end for end for Note that Eis a vector of two elements, one for the each axis of the rectangular coordinate system. The E- eld can also be expressed in terms of its magnitude and the associated angle.

Accuracy of numerical solution depends on the grid size xand y.

A coarse grid (i.e. large xand y) will lead to fast computation of the numerical solution, however will not be accurate. On the other han d a ner grid will require long time for computation but can result in a more accu rate solution. If analytical solution is not known, one can solve the proble m for using successively and compare solutions. An acceptable solution is o btained when the solution does not change signi cantly if grid size is reduced a ny further.

You can test accuracy of your Matlab code for the following problem that can be solved analytically: Write a Matlab code to nd the electric potential V(x, y ) and the electric eld E(x, y ) at all grid points in the domain.

Compare your numerical solution and the analytical solution (given b elow).

Consider the rectangular domain as shown with the applied voltage on the boundary: 3 Solution7 x y V=10 Volts V=0 V=0 V=0 2 1 Fig 3. Domain Analytical solution: It can be shown that the analytical solution of this problem is given by V(x, y ) = m =1 ,3 ,5 ,··· 4 V 0 m 1 sinh m a b sinh m x b sin m y b (14) where a= 2 is the width and b= 1 is the height of the region, and V 0 = 10 is the applied voltage. 4 Project8 4 Project Consider the rectangular domain as shown with the applied voltage onthe boundary: x y V=10 Volts V=0 V=0 2 1 V = 10 Volts Fig 4. Domain • Write a Matlab code to nd the electric potential V(x, y ) and the electric eld E(x, y ) at all grid points in the domain.

• Draw equipotential lines in the domain.

• Draw contour plot of constant eld intensity.

• How much energy is stored in this device?

• Could this device be considered as a capacitor? If so, nd its capacitance.

• Disucss anything else about this problem that you nd interesting or important.