We proceed now with the solution of Equation 50 on the basis of the Finite Element Method [29] and [30]. As a first step we introduce the residual

which must be zero in accordance with the state problem. However, it is impractical to enforce *R*( *x *) = 0 at every point in the domain from *x *= 0 to *x *= *x*_{f }. Since φ( *x*_{f}* *) is not expected to vary substantially over a small distance, say Δ*x *, we subdivide the domain into small segments and instead enforce the condition

over each of the segments. Remarking that *W*_{m}* *( *x *) is some weighting function to be defined later, Equation 84 enforces the differential equation on an average sense over the *N*_{e} *th *segment. By changing the integration or testing interval (and weighting function *W*_{m }) from *m *= 1 to *m *= *N*_{e }, we can construct *N** _{e} *equations for the solution of the discretized field values. Before proceeding to do so, we make the following observations about

*W*

_{m}

*(*

*x*) :

• If *W**m** *( *x *) = δ( *x *– *x*_{m}* *) or *W*_{m}* *( *x *) = δ[ *x *– ( *x*_{m }+ 1 + *x*_{m }) ⁄ 2 ] the resulting weighted residual procedure is referred to as point matching and leads to a form of the finite difference method.

• If *W*_{m}* *( *x *) is set equal to the basis functions used for the representation of φ( *x *) , the procedure is referred to as Galerkin’s method. This is the most popular testing/weighting method for casting the differential equation to a linear system.

• The choice of *W*_{m}* *( *x *) is not completely arbitrary. For the mathematical steps in the FEM procedure to hold rigorously, *W*_{m}* *( *x *) and its derivative must be at least square integrable over the domain. Specifically, for the problem at hand, it must satisfy the condition

In addition, *W*_{m}* *( *x *) must satisfy conditions at the boundary nodes (end-points at *x *= 0 and *x *= *x*_{m }) for the one-dimensional problem which are compatible with the imposed boundary conditions. Certain smoothness conditions on *W*_{m}* *( *x *) and may also need to be imposed.

Before generating a linear of equations from Equation 84 subject to the boundary conditions, it is first necessary to cast in a more suitable form by following the steps:

Step 1: Take advantage of the weighting *W*_{m}* *( *x *) to reduce the order of the derivatives contained in *R*( *x *) . To do so we employ integration by parts, giving:

The last right hand side term can be evaluated by enforcing the known boundary conditions at the endpoints. Its effect on the overall system will be considered later.

Step 2: Derive the “weak form” of the differential equation. The weak form of the differential equation is most appropriate for numerical solution and is obtained by substituting Equation 86 into Equation 84.

which holds providing Equation 50 is valid. However, because of the integral, the weak form Equation 87 enforces the differential equation on an average (and therefore weaker) sense.

**a) Discretization of the “weak” differential equation**

The discretization of Equation 87 to a linear set of equation is done by introducing an expansion for φ( *x *) and then making appropriate choices for the weighting functions. We chose the linear representation

where φ^{e}_{j }( *x *) are the unknown coefficients of the expansion. When *e**o *= 2 , we define linear elements. In this case, the basis functions for element *e *are defined as

**Figure 5: Nodal expansion function for e ^{th} functions considering linear approximation**

The basis functions have unit magnitude at one node and vanish at all others with linear variation between the nodes. When *e**o *= 3 , we have quadratic elements that are also known as second order elements. Each element has three nodes, one of two endpoints, and the third is usually placed at the center of the element. Within each element, the unknown function is approximated as a quadratic function

Enforcing Equation 91 at the three nodes of the element yields:

Solving for *a*^{e}*, b*^{e}*, and c*^{e }and substituting them back into Equation 91, we obtain

where the interpolation or expansion functions are given by

** Figure 6: The nodal expansion function for e ^{th} functions considering quadratic approximation**

When the expansion Equation 92 is substituted into Equation 87 we get

The terms in the brackets are due to contributions from endpoints of the domain and their evaluation is subject to the specific boundary conditions. This equation now explicitly shows how the boundary conditions enter into the construction of the linear system. Hereon, we will refer to their contributions as endpoints since we have not yet specified the type of boundary condition to be imposed.

We are now ready to make different choices for the weighting function to generate a system of linear equations for the solution of { φ _{n}} . As stated earlier, this step is also referred to as testing and Galerkin’s method is usually employed in the finite element method. Specifically, we choose *W*_{m}* *( *x *) = *N*_{j}^{e }( *x *) and for each of these testing or weighting functions a single linear equation is generated.

From Equation 96 we have

We can rewrite Equation 97 in matrix form as

where { φ^{e}* *} is the element electric or magnetic field, { 0 } is a null vector and the element matrices [ *K*^{e}* *]_{ij} and [ *M*^{e}* *]_{ij }are given by

Here [ *K*^{e}* *]_{ij} and [ *M*^{e}* *]_{ij} are tridiagonal when we use linear base function to expand the field element or pentadiagonal when we adopt quadratic base function to expand the field. Since *N _{i}^{e}*

*(*

*x*) are linear functions, the evaluation of integrals can be carried out in closed form provided

*p*

^{e}

*(x*

*)*and

*q*

^{e}

*(x*

*)*are taken as constant over integration of the

*eth*element. Specifically, setting

*p*

^{e}

*(*

*x*) ≅

*p*

*and*

^{e}*q*

^{e}

*(*

*x*) ≅

*q*

*for*

^{e}*x*

_{l}^{e }<

*x*<

*x*

_{2}^{e }.

As a result of integral evaluation we have

(1) Linear elements:

(2) Quadratic elements:

Here l_{e = }x_{2}^{e }– x_{1}^{e}.

**b****) Boundary Condition**

The endpoint contributions appear only when *e *= 1 or *e *= *m *and vanish when the Neumann

or Dirichlet conditions φ|_{x = 0} = φ| _{x = xf }are imposed.

To set the TBC as defined in “Transparent Boundary Condition” on page 23, we

define:

Here,

where the real parts of *k*_{1}, * _{i} *and

*k*

_{m},

*must be restricted to be positive to ensure only radiation outflow.*

_{i}For Dirichlet, Neumann, and TBC boundary conditions, we set *s *= 1 . For PML, we introduce the parameter *s *as defined by Equation 45.