Return to Article Details A defect-correction nodal finite element method for time-dependent Maxwell's equations on polygonal domains

A Defect-Correction nodal Finite Element Method for Time-dependent Maxwell’s equations on Polygonal DomainsThanks: Department of Mathematics, University of Buea, Molyko, Buea, Cameroon, e-mail: nkeck.jake@ubuea.cm, https://orcid.org/0000-0003-0492-8247.

Jake Léonard Nkeck
Date: January 06, 2026; accepted: April 15, 2026; published online: May 05, 2026.
Abstract.

This paper develops a nodal finite element Crank-Nicolson method of lines to solve the time-dependent Maxwell’s equations on polygonal domains with re-entrant corners. Nodal Finite Element methods are used to solve Maxwell’s equations with an optimal convergence rate when the domain is convex or has a smooth boundary, but may fail to converge if the domain has a re-entrant corner. The Defect-Correction method presented is based on a decomposition of the solution in terms of Fourier and Bessel’s series, an extraction of the singular function and an approximation of the regular part of the solution. Optimal convergence results are recovered using the method in both the energy norm and the L2-norm.

Key words and phrases: 
Maxwell’s Equations, Crank-Nicolson Method, Re-entrant Corner, Singular function
2005 Mathematics Subject Classification
37A17, 76F02, 76M10, 74S20, 74S25

1. Introduction

Electromagnetic waves are of particular importance in applied physics, engineering, and materials science, for example, in radars, antennas, and the detection of cracks in metals. Electromagnetic phenomena are modeled by a particular system of partial differential equations (PDEs) called Maxwell’s equations or equations of electromagnetism. Maxwell’s equations involve a couple of vector functions (u,b) called electromagnetic field, where u is the electric field and b is the magnetic field.

A general result of numerical methods for PDEs is that the accuracy of approximation depends on the regularity of the unknown solution, see [1, 2, 3, 4]. It is well known that solutions of boundary value problems in domains with corners, cracks, edges and conic vertices may entail singularities, even if the data are smooth. Moreover, the asymptotic behavior of the solution near these geometric singularities has been well studied and the results are now well known, see, for example [5, 6, 7, 8]. However, unlike general elliptic boundary value problems, the solution of Maxwell’s equations in domains with geometric singularities exhibits some regularity properties that make its approximation particularly difficult.

In the case of time-harmonic Maxwell’s equations in a two-dimensional domain Ω with corners and for a given right hand side function in the Hilbert space L2(Ω)2, it has been shown that the solution belongs to the Sobolev space Hm(Ω)2 (m=1,2) only if the value of the largest angle is smaller than π/m, see, for example [9, 10, 11, 12]. It follows that the solution belongs to the space H1(Ω)2 only if Ω is convex and to the space H2(Ω)2 only if the largest angle at the corners is less than π/2. We would say that the solution of the boundary value problem has geometric singularities if for a given right hand side function in L2(Ω)2 the solution belongs only to the Sobolev space Hs(Ω)2 with s<2. It follows that on non-convex domains, the standard H1-conforming FEM cannot be employed for the solution of Maxwell’s equations.

The work of Costabel [13] in 2002 combined with the bilinear form introduced in [14] has allowed the scientific community to renew the nodal FEM for Maxwell’s Equations and some new finite element schemes have been developed in order to recover the convergence rate on non-convex domains. We mention in particular:

  1. -

    The singular field method (SFM) introduced by Hazard and Lenoir [15], Hazard et al. [16, 17]. This method consists in splitting the solution u into two parts u=ur+γsus, with a ”regular part” ur that lies in the usual Sobolev space H1 and a singular part us that lies in a known finite dimensional space and a constant γs to be determined from the finite element linear system. The constant γs is known as the coefficient of singularity. With this method one can approximate the solution in non-convex polygonal domains with nodal FEM but the rate of convergence is still lower than the one on smooth domains because the splitting is not optimal.

  2. -

    The Orthogonal Singular Field Method (OSFM) introduced by Assous et al. [17], where the solution u is again split into two parts u=ur+(ur+γsus), a regular part ur that lies in a subspace of H1 and the second part (ur+us) is in a subspace of the space of solutions orthogonal to the regular space. The OSFM and the SFM have the same rate of convergence but the OSFM is more stable due to the omission of the use of cut-off functions present in the SFM for the instabilities due to the use of cut-off functions.

  3. -

    The ”λ-approach” for two-dimensional vector problems developed by Jamelot in 2004 [18], where the splitting u=ur+γsus is done but γs is given by a formula and depends only on the domain and the initial data. This method ameliorates the rate of convergence of the SFM and the OSFM but the rate is still not optimal (of order 𝒪(h2π/ω1ε) in the energy norm and 𝒪(h4π/ω2ε) in the L2-norm, ω being the greatest value of the angles of the domain).

  4. -

    The weighted regularization developed by Costabel and Dauge [13], where the boundary conditions are multiplied by a weight that depends on the distance to the geometric singularities. In this method each singularity requires the construction of a new weight. It is stated in [19] that this method requires the approximation space to contain the gradient of C1-scalar functions, which excludes low order finite element spaces (C0-finite element spaces for example). This restriction is removed in the work of Buffa et al. [20] by considering a mixed form of the weighted L2-stabilization technique on special meshes. This method has also been simplified by Otin in [21, 22] where the method is performed by using a weight equal to zero in the elements near the singularity and equal to one in the other elements.

  5. -

    The mixed methods with natural boundary conditions developed by Ciarlet Jr. et al.[4]. These methods consist in the dualisation of the equation on the divergence and the relation on the tangential or normal trace of the field with some Lagrange multipliers. Several authors have also contributed to the development of these methods. We highlight the works of Codina et al. [23, 24, 25] where a novel augmented formulation is produced by adding the Laplacian of the Lagrange multiplier multiplied by a mesh dependent stabilizing term to the equation resulting from the dualisation of the divergence equation.

  6. -

    The L2-projection methods developed by Duan and al. [26, 27]. In these methods, L2-projectors are applied to both curl and div formulations and linear continuous finite elements enriched with some higher order bubble functions are employed in order to approximate low regular functions. These methods do not impose information on the geometric singularities of the domain boundary but are still limited to linear continuous elements.

  7. -

    The interior penalty method [19, 28], where the idea consists of controlling the divergence of the electric field in a Sobolev space with fractional negative exponent. The optimal rate of convergence is recovered with this method.

  8. -

    The Predictor-Corrector nodal FEM that makes use of the explicit extraction formulas for the coefficients of the singularities of the solution near the corners ( see [9, 11, 12, 29, 30]). The optimal convergence rate is recovered by the method. The method can be applied with high order finite element polynomials but faces the problem of logarithmic singularities that may easily happen in the case of Maxwell’s equations.

This paper extends the Predictor-Corrector nodal FEM developed for the 2D time-harmonic Maxwell’s equations on convex and non-convex polygonal domains with exactly one re-entrant angle Ω centered at the origin O(0,0). The method is based on a Fourier decomposition of the solution and a defect-correction algorithm to recover the optimal convergence known for Maxwell’s equations on convex polygonal domains or on domains with a 𝒞2-boundary.

This paper is organized as follows, Section 2 presents the Maxwell’s equations on polygonal domains, Section 3 develops a local decomposition of the solution around the singular corner, Section 4 proposes a Defect-Correction FEM with the error estimates in L2- and the energy norm and Section 5 concludes the paper.

2. The Model Problem

Given a vector function v=(v1,v2)T, define the divergence of v by

div v= v1x1+v2x2.
curlv= v2x1v1x2.

Given a scalar function v, define curl v by

curl v= (vx2,vx1).

Define L2(Ω) to be the space of classes of measurable and square integrable functions over Ω equipped with the norm u0=(Ω|u(x1,x2)|2dx)1/2 that will also denote the norm in the Cartesian product space L2(Ω)2, dx=dx1dx2 being the Lebesgue measure.

H0(curl,Ω) :={vL2(Ω)2:curlvL2(Ω) and vn=0 on Ω}
H(div,Ω) :={vL2(Ω)2:div vL2(Ω)}
H(div0,Ω) :={vL2(Ω)2:div v=0 in Ω}
H0(curl,div,Ω) :=H0(curl,Ω)H(div,Ω)
Hm(Ω) :={v:α1+α2vx1α1x2α2L2(Ω),α1,α2,α1+α2m}
HN(Ω) :={vH1(Ω)2:vn=0 on Γ}
H01(Ω) :={uH1(Ω):u=0 on Γ}

equipped with the norms

vcurl :=(v02+curlv02)1/2
vdiv :=(v02+div v02)1/2
vcd :=(v02+curlv02+div v02)1/2
vm :=(v02+α1,α2mα1+α2vx1α1x2α202)1/2.
vHN(Ω) :=v1=:vH01(Ω)

where m denotes the norm in Hm(Ω) and also the norm in the cartesian product space Hm(Ω)2. The quantity defined by

|v|cd:=(curlv02+div v02)1/2

is the semi-norm in H0(curl,div,Ω).

Given an interval I and a Hilbert space X, 𝒞k(I,X) will denote the space of bounded k-continuously differentiable functions u on I of the form tu(,t)X equipped with the norm

u𝒞k=suptI,|α|kαu(,t)X.

The space Ck(I,X) will denote the space of functions u:=(u1,u2), u1,u2𝒞k(I,X) equipped with the norm

uCk=(u1𝒞k2+u2𝒞k2)1/2.

The space L2(I,X) is the completion of 𝒞0(I,X) with respect to the norm

uL2(I,X)=(Iu(,t)X2𝑑t)12.

The space L2(I,X) is the space of functions u:=(u1,u2), u1,u2L2(I,X) equipped with the norm

uL2(I,X)=(u1L2(I,X)2+u2L2(I,X)2)1/2.

The space Hm(I,X), m, m>0 is the space of functions tu(,t)X such that α1+α2utα1+α2L2(I,Ω) with the norm

uHm(I,X)=(Iα1+α2mα1+α2utα1+α2X2𝑑t)12.

The space Hm(I,X) is the space of functions u:=(u1,u2), u1,u2Hm(I,X) equipped with the norm

uHm(I,X)=(u1Hm(I,X)2+u2Hm(I,X)2)1/2.

Given a simply connected polygonal domain Ω2 with boundary Ω, a function f:=(f1,f2)TH1([0,T],L2(Ω)2) such that div f=0, two functions u0:=(u10,u20)TH0(curl,div,Ω)H(div0,Ω)H2(Ω)2 and u1:=(u11,u21)TH(div0,Ω), find u:=(u1,u2)TL2(0,T,H0(curl,div,Ω)) such that

(1) κ22ut2+curlcurlu= f,inΩ×(0,T),
divu= 0,inΩ×(0,T),
un= 0,onΩ×(0,T),
u(,0)= u0(),inΩ,
ut(,0),= u1(),inΩ,

where n=(n1,n2) is the unit outward normal to Ω, T>0 and κ.

Using the formula curl curlu=Δu+(div u) one obtains the system

κ22ut2Δu= f,inΩ×(0,T),
divu= 0,inΩ×(0,T),
(2) un= 0,onΩ×(0,T),
u(,0)= u0(),inΩ,
ut(,0)= u1(),inΩ,

Problem (2) is equivalent to an initial boundary value problem useful for the Lagrange nodal finite element investigation due to the continuity of the finite element functions along the tangential and the normal components. That result is given by the following proposition.

Proposition .

If κ2 is not an eigenvalue of the Laplace operator Δ(), Problem (2) is equivalent to the following problem:

find uL2(0,T,H0(curl,div,Ω)) such that

κ22ut2Δu= f,inΩ×(0,T),
divu= 0,onΩ×(0,T),
(3) un= 0,onΩ×(0,T),
u(,0)= u0(),inΩ,
ut(,0)= u1(),inΩ,
Proof.

A solution of Problem (2) is obviously a solution of Problem (2).

Suppose u is a solution of Problem (2) and let φ=div u, taking the divergence in the first equation of (2), one obtains that φL2(0,T,H01(Ω)) is the unique solution of the Problem

κ22φt2Δφ =0,inΩ×(0,T),
(4) φ =0,onΩ×(0,T),
φ(,0) =div u0()=0,inΩ,
φt(,0) =div u1()=0,inΩ,

But Problem (2) has the unique null solution provided that κ2 is not an eigenvalue of the Laplace operator Δ(), then φ=0 in Ω×(0,T) and u is also a solution of Problem (2). ∎

Taking the inner product of the first, fourth and the fifth equation of Problem (2) with vH0(curl,div,Ω), the integral over Ω and using integration by parts lead to the variational problem: find uL2(0,T,H0(curl,div,Ω)) such that

κ2d2dt2Ωuvdx+a(u,v)= Ωfvdx,vH0(curl,div,Ω)
(5) Ωu(,0)vdx= Ωu0vdx,vH0(curl,div,Ω)
Ωut(,0)vdx= Ωu1vdx,vH0(curl,div,Ω)

where a(u,v):=Ωcurlucurlv+divudivvdx. The following results prove the existence and uniqueness of the solution of the variational problem (2). To prove it, one can apply Theorem 8.1, page 287 of [31].

Proposition \theproposition.

If fH1(0,T,H0(curl,divΩ)), u0H0(curl,div,Ω)H(div0,Ω)H2(Ω)2 and u1H(div0,Ω), then the variational problem (2) has a unique solution uC0([0,T],H0(curl,div,Ω)) such that ut belongs to C0([0,T],H(div0,Ω)). Furthermore the solution depends continuously on the data.

3. Decomposition of the Solution

An approximation of the solution of the variational problem (2) involves a particular study of the relation between HN(Ω) and H0(curl,div,Ω). From [16] page 2033, H0(curl,div,Ω)=HN(Ω)XS where

XS:={φH01(Ω):f𝒩:Ωφψdx=ΩfψdxψH01(Ω)},

𝒩 being the orthogonal of Δ(H2(Ω)H01(Ω)) in L2(Ω).

This introduction allows the decomposition of the solution into a space regular and a space singular part. The geometric singularity being a local problem see [8, p. 71]it is judicious to study the problem in a circular sector.

For sake of simplicity, it is assumed that the domain Ω has only one corner centered at the origin O(0,0) in the (x1,x2)-plane with angle ω>π/2, ωπ. Introduce the polar coordinates (r,θ) with x1=rcosθ, x2=rsinθ, 0θω. Let R0>0, consider the restriction on a circular sector G0 with

G0¯={(rcosθ,rsinθ):0θω,0rR0},

and the cut-off function

η(r):={1,if 0r<R0/30η(r)1,if R0/3r2R0/30,if r>2R0/3.

Set uη:=ηu, uη describes the solution u around the corner O(0,0) and uη is the unique solution of the problem: find uηL2(0,T,H0(curl,div,G0)) such that

(6) κ22uηt2Δuη= fη,inG0×(0,T),
(7) divuη= 0,onG0×(0,T),
(8) uηn= 0,onG0×(0,T),
(9) uη(,0)= ηu0(),inG0,
(10) uηt(,0)= ηu1(),inG0,

where

fη:=(ηf1u1Δη2ηu1ηf2u2Δη2ηu2).

The one-to-one mapping (x1,x2)(r,θ) transforms G0 into a rectangle G~0:={(r,θ):0rR0,0θω} and we can pass through polar coordinates by setting u~(r,θ,t)=uη(x1,x2,t) and f~(r,θ,t)=fη(x1,x2,t)

(uruθ)= (uη1cosθ+uη2sinθuη1sinθ+uη2cosθ),(frfθ)=(fη1cosθ+fη2sinθfη1sinθ+fη2cosθ)
(ur0u0θ)= (ηu10cosθ+ηu20sinθηu10sinθ+ηu20cosθ),(ur1u1θ)=(ηu11cosθ+ηu21sinθηu11sinθ+ηu21cosθ),

Problem (6) becomes

(11) κ22urt22urr21r22urθ21rurr+2r2uθθ+1r2ur= fr,inG0~×(0,T)
κ22uθt22uθr21r22urθθ21ruθr2r2urθ+1r2uθ= fθ,inG0~×(0,T)
urr+1rur+1ruθθ= 0,ur=0,ifθ=0
urr+1rur+1ruθθ= 0,ur=0,ifθ=ω
|ur(0,θ,t)|<,|uθ(0,θ,t)|< ,tT
ur(R0,θ,t)=uθ(R0,θ,t)= 0, 0<θ<ω 0tT
ur(,0)=ur0(),uθ(,0)= uθ0()
urt(,0)=ur1(),uθt(,0)= uθ1()

The derivatives are interpreted in the sense of distributions. The boundary conditions at θ=0 and θ=ω allow us to consider the following Fourier decompositions for ur and uθ,

ur(r,θ,t) =k=1urk(r,t)sinλkθ, uθ(r,θ,t)=k=1uθk(r,t)cosλkθ
fr(r,θ,t) =k=1frk(r,t)sinλkθ, fθ(r,θ,t)=k=1fθk(r,t)cosλkθ
ur0,1(r,θ) =k=1urk0,1(r)sinλkθ, uθ0,1(r,θ)=k=1uθk0,1(r)cosλkθ.

where λk=kπ/ω, k, k>0. The system (11) becomes

(12) κ22urkt22urkr2+λkr2urk1rurkr2λkr2uθk+1r2urk =frk
κ22uθkt22uθkr2+λkr2uθk1ruθkr2λkr2urk+1r2uθk =fθk
urkr+1rurkλkruθk=0,urk =0,ifθ=0 or θ=ω
|urk(0,t)|<,|uθk(0,t)| <,0tT
urk(R0,t)=uθk(R0,t) =0,0tT
urk(r,0)=η(r)urk0(r),uθk(r,0) =ηuθk0(r),0rR0
urkt(r,0)=η(r)urk1(r),uθkt(r,0) =ηuθk1(r),0rR0

Setting u3k=urk+uθk and u4k=urkuθk the first and the second equations of Problem (12) become

(13) κ22u3kt22u3kr2+λk2+1r2u3k1ru3kr2λkr2u3k=frk+fθkκ22u4kt22u4kr2+λk2+1r2u4k1ru4kr2λkr2u4k=frkfθk

The homogeneous equations associated to the equations (13) can be solved by separation of variables by setting uik(r,t)=φik(t)ψik(r), i=3,4 to obtain the equality φik′′φik=1κ2ψik(ψik′′(λk1)2r2ψi+1rψik), i=3,4.

Due to the mixed boundary conditions at t=0 from the two last equalities of (12) and the smoothness of the solution u in time, the operator d2()dt2 has positive and discrete eigenvalues αkm2, m, m1 arranged in an increasing sequence. One writes φik′′φik=1κ2ψik(ψik′′(λk1)2r2ψi+1rψik)=αkm2, i=3,4, m, so, the equations ψik′′+1rψik((λk1)2r2κ2αkm2)ψi=0, i=3,4 lead to ψik(r)=m=1C1imJλk1(καkmr)+C2imYλk1(καkmr), i=3,4 where C1im, C2im are constants and the αkm, m, m1 form an increasing sequence of positive numbers such that Jλk1(καkmR0)=Yλk1(καkmR0)=0, Jλk1 and Yλk1 are the Bessel functions of the first and second kind respectively.

The boundary conditions |urk(0,t)|< and |uθk(0,t)|< in Problem (12) imply that C2im=0, i=3,4, m, m1. Then

urk(r,t)=m=1urkm(t)Jλk1(καkmr), uθ(r,t)=m=1uθkm(t)Jλk1(καkmr)
frk(r,t)=m=1frkm(t)Jλk1(καkmr), fθk(r,t)=m=1fθkm(t)Jλk1(καkmr)
urk0(r)=m=1urkm0Jλk1(καkmr), uθk0(r)=m=1uθkm0Jλk1(καkmr)
urk1(r)=m=1urkm1Jλk1(καkmr), uθk1(r)=m=1uθkm1Jλk1(καkmr)

where

frkm(t)= 1Jλk1(καkm)020R0frk(r,t)Jλk1(καkmr)r𝑑r
= 2ωJλk1(καkm)020ω0R0fr(r,θ,t)Jλk1(καkmr)sin(λkθ)r𝑑r𝑑θ,
fθkm(t)= 1Jλk1(καkm)020R0fθk(r,t)Jλk1(καkmr)r𝑑r
= 2ωJλk1(καkm)020ω0R0fθ(r,θ,t)Jλk1(καkmr)cos(λkθ)r𝑑r𝑑θ,
urkm0= 1Jλk1(καkm)020R0η(r)urk0(r)Jλk1(καkmr)r𝑑r
= 2ωJλk1(καkm)020ω0R0η(r)ur0(r,θ)Jλk1(καkmr)sin(λkθ)r𝑑r𝑑θ,
urkm1= 1Jλk1(καkm)020R0η(r)urk1(r)Jλk1(καkmr)r𝑑r
= 2ωJλk1(καkm)020ω0R0η(r)ur1(r,θ)Jλk1(καkmr)sin(λkθ)r𝑑r𝑑θ.

Replacing urk and uθk by their decompositions in Problem (12) leads to

(14) urkm(t)=(1αkmurkm1+1κ2αkm0tfrkm(τ)cos(αkmτ)dτ)sin(αkmt)+(urkm01κ2αkm0tfrkm(τ)sin(αkmτ)dτ)cos(αkmt),uθkm(t)=(1αkmuθkm1+1κ2αkm0tfθkm(τ)cos(αkmτ)dτ)sin(αkmt)+(urkm01κ2αkm0tfθkm(τ)sin(αkmτ)dτ)cos(αkmt).

Then

ur(r,θ,t)=k,m=1urkm(t)Jλk1(καkmr)sinλkθ,uθ(r,θ,t)=k,m=1uθkm(t)Jλk1(καkmr)cosλkθ,

where urkm and uθkm are defined in (14). Using the fact that

Jλk1(καkmr)=n=0(1)nn!Γ(n+λk)(καkmr2)2n+λk1,

the solution u~:=(ur,uθ)T has the decomposition

(15) u~(r,θ,t) = k,m=1n=0sk,m,n(r,θ,t)
= w~(r,θ,t)+k,m=1n=0,2n+λk1<2s~k,m,n(r,θ,t)
= w~M(r,θ,t)+k=1n=0,2n+λk1<2m=1Ms~k,m,n(r,θ,t),

where M, M1,

s~k,m,n(r,θ,t)=(1)nn!Γ(n+λk)(καkmr2)2n+λk1(urkm(t)sinλkθuθkm(t)cosλkθ),

and wC1(0,T,H2(Ω)2).

Remark .

In cartesian coordinates, uη has the decomposition

(16) uη(x1,x2,t):=w(x1,x2,t)+k=1n=0,2n+λk1<2m=1Msk,m,n(r,θ,t),

where

sk,m,n(x1,x2,t)=(1)nn!Γ(n+λk)(καkmr2)2n+λk1(urkm(t)sin(λk1)θuθkm(t)cos(λk1)θ).
Lemma .
(17) |k,m=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1urkm(t)|<

and

(18) |k,m=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1uθkm(t)|<
Proof.
|k,m=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1urkm(t)|
k,m=1n=0,2n+λk1<2|καkm|24n!|Γ(n+λk)||urkm0cos(αkmt)
+1καkmurkm1sin(αkmt)+1κ2αkm0tfrkm(τ)sin(αkmτ)dτ|.

The fact that rJλk1(καkmr)=1καkm(λkJλk(καkmr)+rdJλkdr(καkmr)) implies that

0R0ur0(r)sin(λkθ)Jλk1(καkmr)𝑑r=
= λkκαkm0R0ur0(r)sin(λkθ)Jλk(καkmr)dr
+1καkm0R0ur0(r)sin(λkθ)dJλkdr(καkmr)rdr
= λkκαkm0R0ur0(r)sin(λkθ)rλk1(r1λkJλkdr(καkmr))dr
+[rκαkmur0(r)sin(λkθ)Jλk(καkmr)]0R0
1καkm0R0(ur0(r)+rdur0dr(r))sin(λkθ)Jλk(καkmr)dr
= λkκαkm0R0ur0(r)sin(λkθ)rλk1ddr(r1λkJλk1(καkmr))dr
+1κ2αkm20R0(ur0(r)+rdur0dr)sin(λkθ)rλk1(καkmr1λkJλk(καkmr))dr
= 1καkmλkκ2αkm0R0ur0(r)sin(λkθ)rλk1ddr(r1λkJλk1(καkmr))𝑑r
+1κ2αkm20R0dur0drsin(λkθ)rλkddr(r1λkJλk1(καkmr))dr
= 1κ2αkm2[((1καkm)ur0+rdur0dr)Jλk1(καkmr)]0R0
1κ2αkm20R0((καkmλk1)2rur0+dur0dr+rd2ur0dr2)sin(λkθ)Jλk1(καkm)dr
= 1κ2αkm20R0((καkmλk1)2r2ur0+1rdur0dr+d2ur0dr2)sin(λkθ)Jλk1(καkm)rdr
CJλk1(καkm)0|καkm|2u02, if u0H2(G0)2.

One also has

0R0ur1(r,θ)sin(λkθ)Jλk1(καkmr)r𝑑r=
= 1καkm0R0ur1(r)sin(λkθ)Jλk1(καkmr)καkmr𝑑r
= 1καkm0καkmR0ur1(Rκαkm,θ)sin(λkθ)Jλk1(R)R𝑑Rwhere R=καkmr
1|καkm|(0καkmR0|ur1(Rκαkm,θ)R1/2|2𝑑R0καkmR0|Jλk1(R)R1/2|2𝑑R)12
CJλk1(καkm)0u10.

Using the same way one can show that

0R0frkm(r)sin(λkθ)Jλk1(καkmr)drCJλk1(καkm)0f0.

Hence (3) becomes

|k,m=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1urkm(t)|
Ck,m=1n=0,2n+λk1<2[(u02+1καkmu10)2ωn!|Γ(n+λk)|Jλk1(καkm)0
+1|καkm|20tsin(αkmτ)dτf0]
|k,m=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1urkm(t)|
k,m=1n=0,2n+λk1<2C(u02+1|καkm|u10+Tf0|καkm|2)2ωn!|Γ(n+λk)|Jλk1(καkm)0
k,m=1n=0,2n+λk1<2C(R0,T,u0,u1,f)2ωn!|Γ(n+λk)|Jλk1(καkm)0
< .

Similar methods can be used with uθ0, uθ1 and fθkm to show (18). ∎

A truncation of the summation on the integer m can be done and the following error estimate is obtained.

Lemma .

Let M, M>1, and let {αkm}m1 be an increasing sequence of positive numbers such that Jλk1(καkmR0)=0 where λk=kπ/ω, R0>0 is fixed and κ is defined as in (1). Let

cr(t)=k,m=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1urkm(t),
cθ(t)=k,m=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1uθkm(t),
cr,M(t)=m=1Mk=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1urkm(t),
cθ,M(t)=m=1Mk=1n=0,2n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1uθkm(t).

Then

(20) |cr(t)cr,M(t)|Cmaxk1,n0, 2n+λk1<2{|αkM|λk}

and

(21) |cθ(t)cθ,M(t)|Cmaxk1,n0, 2n+λk1<2{|αkM|λk}.
Proof.
|cr(t)cr,M(t)|= |m=M+1k=12n+λk1<2(1)nn!Γ(n+λk)(καkm2)2n+λk1urkm(t)|
m=M+1k=12n+λk1<2C2ωn!|Γ(n+λk)|Jλk1(καkm)0

But

Jλk1(καkmr)=(καkmr)λk12λk1Γ(λk)(1(καkmr)22(2λk)+(καkmr)42(4)(2λk)(2λk+2))

then

|cr(t)cr,M(t)| C(R0)m=M+1k=1n=0,2n+λk1<22λk1Γ(λk)2ωn!|Γ(λk+n)||αkm|λk
Cmaxk1,n0, 2n+λk1<2{|αkM|λk}

The inequality (21) can be deduced using the same manner. ∎

Remark .

αkM as M but since λk>0, |αkM|λk0 as M.

4. The Defect-Correction Finite Element method

4.1. A Semi-discrete Approximation

This section presents a space discretization of the solution by a Galerkin approximation with H1-conforming finite element methods. The Galerkin method considers a quasi-uniform and shape-regular triangulation 𝒯h of Ω with diameter h, and finds a semi-discrete approximation uhL2(0,T,Vh) of the solution uL2(0,T,H(curl,div,Ω)) of Problem (2) satisfying

(22) κ2d2dt2Ωuhvhdx+a(uh,vh)=Ωfvhdx,vhVhΩuh(,0)vhdx=Ωu0vhdx,vhVhΩuht(,0)vhdx=Ωu1vhdx,vhVh

where Vh:={vh𝒞(Ω)2:vh|K𝒫1(K)2and vhn=0 on Ω}HN(Ω) is a finite dimensional space.

The following result presents the standard error estimates for regular solutions uL2(0,T,H2(Ω)2) which happens when Ω is convex or has a 𝒞2-boundary. A proof can be adapted from the one in [32, Theorem 4.1, p.  575],

Theorem 1.

If Ω is a convex polygonal domain with all its angles less than or equal to π/2, u(,t) the solution of Problem (2) and uh(,t) the solution of Problem (22), t(0,T), if uL(0,T,H(curl, div,Ω)) and utL2(0,T,H(curl, div,Ω)) and kutkL2(0,T,L2(Ω)2), k=3,4, then

u(,t)uh(,t)cdChu(,t)uh(,t)0Ch2.

When the domain Ω is convex with an angle greater than π/2 the solution may converge slowly and if Ω is not convex, the solution may fail to converge. The following algorithm based on the decomposition (16) proposes a way to recover the optimal convergence as in Theorem 1.

Algorithm 1 Defect-Correction Algorithm
  1. 1)

    Solve Problem (22) in Vh and obtain the solution uh(0).

  2. 2)

    Compute the sk,m,n(0) using formula (16), replacing u with uh(0), k,m,n, k1, 1mM, M, 2n+λk1<2, that i,s

    sk,m,n(0)(x1,x2,t)=(1)nn!Γ(n+λk)(καkmr2)2n+λk1(urkm(0)(t)sin(λk1)θuθkm(0)(t)cos(λk1)θ),

    where

    urkm(0)(t)= (1αkmurkm(0)1+1κ2αkm0tfrkm(0)(τ)cos(αkmτ)𝑑τ)sin(αkmt),
    +(urkm(0)01κ2αkm0tfrkm(0)(τ)sin(αkmτ)𝑑τ)cos(αkmt)
    uθkm(0)(t)= (1αkmuθkm(0)1+1κ2αkm0tf(0)θkm(τ)cos(αkmτ)𝑑τ)sin(αkmt)
    +(urkm(0)01κ2αkm0tf(0)θkm(τ)sin(αkmτ)𝑑τ)cos(αkmt),
    frkm(0)(t)= 2ωJλk1(καkm)020ω0R0fr(0)(r,θ,t)Jλk1(καkmr)sin(λkθ)r𝑑r𝑑θ,
    fθkm(0)(t)= 2ωJλk1(καkm)020ω0R0fθ(0)(r,θ,t)Jλk1(καkmr)cos(λkθ)r𝑑r𝑑θ,
    (fr(0)fθ(0))= (fη1(0)cosθ+fη2(0)sinθfη1(0)sinθ+fη2(0)cosθ),
    fη(0):= (fη1(0)fη2(0))=(ηf1u(0)1Δη2ηu(0)1ηf2u(0)2Δη2ηu(0)2),
    (ur(0)uθ(0))= (uη1(0)cosθ+uη2(0)sinθuη1(0)sinθ+uη2(0)cosθ),
    urkm(0)0= 2ωJλk1(καkm)020ω0R0η(r)ur0(r,θ)Jλk1(καkmr)sin(λkθ)r𝑑r𝑑θ,
    urkm(0)1= 2ωJλk1(καkm)020ω0R0η(r)ur1(r,θ)Jλk1(καkmr)sin(λkθ)r𝑑r𝑑θ,
    uθkm(0)0= 2ωJλk1(καkm)020ω0R0η(r)uθ0(r,θ)Jλk1(καkmr)cos(λkθ)r𝑑r𝑑θ,
    uθkm(0)1= 2ωJλk1(καkm)020ω0R0η(r)uθ1(r,θ)Jλk1(καkmr)cos(λkθ)r𝑑r𝑑θ.
  3. 3)

    Find whM(1) solution of

    (23) κ2d2dt2ΩwhM(1)vhdx+a(whM(1),vh)=
    =Ω(fκ2k,n,m=1M2sk,m,n(0)t2)vhdxk,m,na(sk,m,n(0),vh),vhVh,
    ΩwhM(1)(,0)vhdx=Ω(u0k,n,m=1Msk,m,n(0)(,0))vhdx,vhVh,
    ΩwhM(1)t(,0)vhdx=Ω(u1k,n,m=1Msk,m,n(0)t(,0))vhdx,vhVh,

    where a(u,v)=Ωcurlucurlv+divudivvdx.

  4. 4)

    Compute uhM(1)=whM(1)+k,n,m=1Msk,m,n(0).

Error Estimates

From the decomposition (16) and the Defect-correction Algorithm 1, one can derive the following error estimates.

Lemma .

Let 𝒯h be a quasi-uniform and shape-regular triangulation of Ω and let Πh be the usual Lagrange interpolation operator over the triangulation 𝒯h. Then for k,m,n, k1, m1, 2n+λk1<2, λk=kπ/ω,

(24) sk,m,n(,t)Πhsk,m,n(,t)cd Chλk1,
(25) sk,m,n(,t)Πhsk,m,n(,t)0 Chλk.
Proof.

Let us divide the triangulation 𝒯h of the domain Ω into two parts h1={K𝒯h:(0,0)K} and h2={K𝒯h:(0,0)K}, one has

sk,m,n(,t)Πhsk,m,n(,t)02 K𝒯hsk,m,n(,t)Πhsk,m,n(,t)02
Kh1sk,m,n(,t)Πhsk,m,n(,t)2
+Kh2sk,m,n(,t)Πhsk,m,n(,t)02.

By the Cauchy-Schwartz inequality

Kh1sk,m,n(,t)Πhsk,m,n(,t)02CKh1sk,m,n(,t)02+Πhsk,m,n(,t)02.

2n+λk1<2 happens when n=0,1, then for Kh1,

sk,m,n(,t)02 C0hK(r2λk1+r2λk+1)(|urkm(t)|2+|uθkm(t)|2)𝑑r
Cmax{hK2λk,hK2λk+2}
Ch2λk.

Πhsk,m,n(,t) being a polynomial of degree 1 then for Kh1,

Πhsk,m,n(,t)02 Cmax{|sk,m,n(x1,x2,t)|2,(x1,x2)K}meas(K)
Ch2λk2hK2
Ch2λk.

Hence Kh1sk,m,n(,t)Πhsk,m,n(,t)02Chλk.

For Kh2, sk,m,n(,t)H2(K)2 then

Kh2sk,m,n(,t)Πhsk,m,n(,t)02Ch4.

Then sk,m,n(,t)Πhsk,m,n(,t)02Cmax{h2λk,h4}Ch2λk and

sk,m,n(,t)Πhsk,m,n(,t)0Chλk.

Using the same method as above, one can prove that

sk,m,n(,t)Πhsk,m,n(,t)cdChλk1.

Lemma .

Let 𝒯h be a quasi-uniform and shape-regular triangulation of Ω and let Πh be the usual Lagrange interpolation operator over the triangulation 𝒯h. Then for M,n,k, m1, k1, n0, 2n+λk1<2, λk=kπ/ω, there exist C>0 such that

(26) m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)0Chλk|αkM|λk

and

(27) m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)cdChλk1|αkM|λk.
Proof.

Let us divide the triangulation 𝒯h of the domain Ω into two parts h1={K𝒯h:(0,0)K} and h2={K𝒯h:(0,0)K}.

m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)02
Kh1m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)2
+Kh2m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)02.

For Kh1, one can use the inequalities (20), (21) and (24) to obtain

m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)0
m=1M(1)nn!Γ(n+λk)(καkmr2)2n+λk1((urkm(t)urkm(0)(t))sinλkθ(uθkm(t)uθkm(0)(t))cosλkθ)0
+ChKλk|αkM|λk
C(hKλkmaxm=1,,M|αkm|λk+hKλk|αkM|λk)
Chλk|αkM|λk.

For Kh2,

m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)0Ch2|αkM|λk

Then

m=1M(sk,m,n(,t)Πhsk,m,n(,t))+m=M+1sk,m,n(,t)02Chλk|αkM|λk

This shows (26).

The inequality (27) can be shown using the same way as previous.

Theorem 2.

Let 𝒯h be a quasi-uniform and shape-regular triangulation of Ω, let M, M1, let uL2(0,T,H0(curl,div,Ω)) be the solution of Problem (2), and uhM(1) be the solution of Problem (22) using the Defect-Correction Algorithm 1, then for any t(0,T),

(28) u(,t)uhM(1)(,t)cd Cmax2n+λk1<2{h,hλk1|αkM|λk},
(29) u(,t)uhM(1)(,t)0 Cmax2n+λk1<2{h2,hλk1|αkM|λk}.
Proof.
uuhM(1) = wk,n,m=1sk,m,nwhM(1)k,n,m=1Msk,m,n(0)
= wwhM(1)+k,n,m=1M[sk,m,nΠhsk,m,n+(Πhsk,m,nsk,m,n(0))]
+k,n,m=M+1sk,m,n

then

uuhM(1)cd
wwhM(1)cd+k,n,m=1Msk,m,nΠhsk,m,ncd+Πhsk,m,nsk,m,n(0)cd
+k,n,m=M+1sk,m,nΠhsk,m,ncd+Πhsk,m,nsk,m,n(0)cd
wwhM(1)cd+Chλk1[(0tuuh(0)02𝑑τ)1/2+|αkM|λk]
C(wwhM(1)cd+hλk1uuhM(0)0+hλk1|αkM|λk)
C(wwhM(1)cd+hλk1wwhM(0)cdCLOSE
+hλk1m=1Msk,m,nΠhsk,m,n0+hλk1m=M+1sk,m,n0
OPEN+hλk1|αkM|λk)
(30) C(wwhM(1)cd+hλk1wwhM(0)cd+hλk1(2+hλk)|αkM|λk)
C(wwhM(1)cd+hλk1wwhM(0)cd+hλk1|αkM|λk).

The next step consists to find an appropriate bound of wwhM(1)cd.

Set Πhw=wh, s=k,m,nsk,m,n and s(0)=k,n,m=1Msk,m,n(0), one has

(31) κ2d2dt2Ωwhvhdx+a(wh,vh)=Ω(fκ22st2)vhdxa(s,vh)vhVh.

Taking the difference of (31) with the first equation of (22) leads to

κ2d2dt2Ω(whwhM(1))vhdx+a(whwhM(1),vh)=
=Ωκ22(ss(0))t2vhdxa(ss(0),vh).

for any vhVh. Then for vh()=2(whwhM(1))t2(,t)Vh, integration by parts gives

12ddtκ2(whwhM(1))t02+12ddt(a(whwhM(1),whwhM(1)))==Ωκ2(ss(0))t(whwhM(1))tdxa(ss(0),(whwhM(1))t).

Taking the integral between 0 and t and the Inequality (27) give

κ2(whwhM(1))t02+a(whwhM(1),whwhM(1))=
= 2[0tΩκ2(ss(0))t(whwhM(1))tdx𝑑τ0ta(ss(0),(whwhM(1))t)𝑑τ]
= 2[0tΩκ2(ss(0))t(whwhM(1))tdxdτa(ss(0),whwhM(1))
+0ta(ss(0),(whwhM(1))t)dτκ22Ω(ss(0))t(ss(0))tdx
12a(ss(0),ss(0))0ta(ss(0),(whwhM(1))t)dτ]
2(ss(0))t0|κ2(whwhM(1))t|0+2ss(0)cd|whwhM(1)|cd
+(ss(0))t0κ2(ss(0))t0+|ss(0)|cd2
C[hλk1|αkM|λk(κ2(whwhM(1))t02+|whwhM(1)|02)+h2(λk1)|αkM|2λk].

Since a(whwhM(1),whwhM(1))=|whwhM(1)|cd2 one deduces that

(κ2(whwhM(1))t02+|whwhM(1)|cd2)1/2 (h2λk2|αkM|2λk1Chλk1|αkM|λk)1/2
Chλk1|αkM|λk.

The semi-norm ||cd being equivalent to the norm cd in H0(curl,div,Ω) (see [18]), one has

(32) whwhM(1)cd(κ2(whwhM(1))t02+|whwhM(1)|cd2)1/2Chλk1|αkM|λk.

Using (32) in (30) leads to

uuhM(1)cd
C[wwhM(1)cd+k,n,m=1M(sk,m,nΠhsk,m,n)+k,n,m=M+1sk,m,ncd]
C(wwhM(1)cd+hλk1|αkMλk|)
(wwhcd+whwhM(1)cd+hλk1|αkMλk|)
C(h+hλk1|αkMλk|)
Cmax2n+λk1<2{h,hλk1|αkM|λk}

Hence, inequality (28) is proved.

Following the same steps used to establish inequality (30) replacing the norm cd by the norm 0 one has

uuhM(1)0
C[wwhM(1)0+k,n,m=1M(sk,m,nΠhsk,m,n)+k,n,m=M+1sk,m,n0]
C(wwhM(1)0+hλk1|αkMλk|)
C(wwh0+whwhM(1)0+hλk1|αkMλk|)
C(h2+hλk1|αkMλk|)
Cmax2n+λk1<2{h2,hλk1|αkM|λk}

and (29) is proved. ∎

Remark .

If M is chosen so that hλk1|αkM|λkh2, the error estimates (28) and (29) are optimal.

4.2. A Crank-Nicolson full discretization

Considering a quasi-uniform and shape-regular triangulation 𝒯h of the polygonal domain Ω and a Galerkin solution uh of the variational problem (22), let {ϕjh}1jN be a basis of the space Vh, one has uh(,t)=j=1Nujh(t)ϕjh()=UhT(t)Φh() where Uh(t)=(ujh(t))1jNT and Φh=(ϕjh)1jN, satisfies the variational system

κ2Mhd2Uhdt2+AhUh =Fh(t),t(0,T),
MhUh(0) =(Ωu0ϕlhdx)1lNT,
MhdUhdt(0) =(Ωu1ϕlhdx)1lNT,

where Mh=(Ωϕjhϕlhdx)1j,lN, Ah=(a(ϕjh,ϕlh))1j,lN and

Fh(t)=(Ωf(,t)ϕlhdx)1lNT.

The full discretization consists to subdivide the interval [0,T] into L subintervals [tj1,tj] of size τ=tjtj1=TL, j=1,,L such that [0,T]=j=1L[tj1,tj], 0=t0<t1<<tL=T. The σ-schemes (0σ1) consists of approximating Uh(tl)=(ujh(tl))1jNT, 0lL such that

(33) κ2τ2Mh(Uh(tl+1)2Uh(tl)+Uh(tl1))+
+σAhUh(tl+1)+(12σ)AhUh(tl)+σAhUh(tl1)
=σFh(tl+1)+(12σ)Fh(tl)+σFh(tl1),
MhUh(0) =(Ωu0ϕlh𝐝𝐱)1lNT,
MhdUhdt(0) =(Ωu1ϕlh𝐝𝐱)1lNT.

So if we set uhτ(,t) be an interpolation of the solution (Uh(tl))1lL, that is, uhτ(,tl)=Uh(tl)Φ() with a Crank-Nicolson scheme (σ=1/2), one has the following error estimates from Theorem 1

u(,tl)uhτ(,tl)cdC(h+τ2),u(,tl)uhτ(,tl)0C(h2+τ2),

for l=1,,L if Ω is a convex polygonal domain and if u belongs to L(0,T,H(curl, div,Ω)), ut is in L2(0,T,H(curl, div,Ω)) and kutk belongs to L2(0,T,L2(Ω)2), k=3,4.

The following algorithm is developed to recover the convergence rate of the Crank-Nicolson scheme when the domain Ω is not convex.

Algorithm 2 Defect-Correction Algorithm
  1. 1)

    Solve Problem (33) and obtain the solution uhτ(0).

  2. 2)

    Compute the sk,m,n(0) using formula (16), replacing u with uhτ(0), k,m,n, k1, 1mM, M, M1, 2n+λk1<2.

  3. 3)

    Find whτM(1) solution of

    (34) κ2d2dt2ΩwhτM(1)vhdx+a(whτM(1),vh)==Ω(fκ2k,m,n2sk,m,n(0)t2)vhdxk,m,na(sk,m,n(0),vh),vhVh,
    Ωw(1)hτM(,0)vhdx=Ω(u0k,m,ns(0)k,m,n(,0))vhdx,vhVh,ΩwhτM(1)t(,0)vhdx=Ω(u1k,m,nsk,m,n(0)t(,0))vhdx,vhVh.
  4. 4

    Compute uhτM(1)=whτM(1)+k,m,nsk,m,n(0).

Error Estimates of the Crank-Nicolson FEM

Theorem 3.

Let 𝒯h be a quasi-uniform and shape-regular triangulation of Ω, let M, M1, let uL2(0,T,H0(curl,div,Ω)) be the solution of Problem (2), and uhτM(1) be the interpolation of the solution of Problem (33) using the Defect-Correction Algorithm 2 in a subdivision of the interval [0,T] in sub-interval [tl1,tl] of size τ=TL, l=1,,L such that [0,T]=l=1L[tl1,tl], 0=t0<t1<<tL=T. Then for any 1lL,

(35) u(,tl)uhτM(1)(,tl)cdC(maxk1,n0, 2n+λk1<2{h,hλk1|αkM|λk}+τ2),
(36) u(,tl)uhτM(1)(,tl)0C(maxk1,n0, 2n+λk1<2{h2,hλk1|αkM|λk}+τ2).
Proof.
u(,tl)uhτM(1)(,tl)cd
u(,tl)uhM(1)(,tl)cd+uhM(1)(,tl)uhτM(1)(,tl)cd
C(max2n+λk1{h,hλk1|αkM|λk})+uh(,tl)uhτM(1)(,tl)cd.

But

uhM(1)(,tl)=whM(1)(,tl)+k,n,m=1Msk,m,n(0)anduhτM(1)(,tl)=whτM(1)(,tl)++k,n,m=1Msk,m,n,τ(0),

​​​​where 𝐬k,m,n,τ(0) is the interpolation of the 𝐬k,m,n(0)(,tl),0lL.

Then

uhM(1)(,tl)uhτM(1)(,tl)cd
whM(1)(,tl)whτM(1)(,tl)cd+k,n,m=1Msk,m,n(0)sk,m,n,τ(0)cd
whM(1)(,tl)whτM(1)(,tl)cd+C(maxk,n, 2n+λk1<2{hλk1|αkM|λk}+τ2).

Also, using (32)

whM(1)(,tl)whτM(1)(,tl)cd
whM(1)(,tl)wh(,tl)cd+wh(,tl)whτM(1)(,tl)cd
whM(1)(,tl)wh(,tl)cd
+wh(,tl)whτ(1)(,tl)cd+whτ(,tl)whτM(1)(,tl)cd
C(maxk,n,2n+λk1<2{hλk1|αkM|λk}+h+τ2+whτ(,tl)whτM(1)(,tl)cd).

But one can use a similar method as the one in Theorem 2 to show that

whτ(,tl)whτM(1)(,tl)cdC(max2n+λk1<2{hλk1|αkM|λk}+τ2).

Hence

u(,tl)uhτM(1)(,tl)cdC(maxk1,n0, 2n+λk1<2{h,hλk1|αkM|λk}+τ2).

One can use the same method as above to show (36). ∎

Remark .

One can remark that if M is chosen such that

maxk1,n0, 2n+λk1<2{hλk1|αkM|λk}h2,

then the optimal convergence of a linear Crank-Nicolson FEM is obtained, that is

u(,tl)uhτM(1)(,tl)cdC(h+τ2)u(,tl)uhτM(1)(,tl)0C(h2+τ2)

5. conclusion

This paper presents a finite element method to solve the time-dependent Maxwell’s equations on non-convex polygonal domains. The nodal FEM is coupled with a Fourier decomposition that extracts the singular behavior of the solution and recovers the optimal linear convergence. The extension to polygonal domains with multiple re-entrant corners can be handled by the method through the use of local cut-off functions near the corners and adding the different singular parts to the final solution. However, the method assumes that the solution is sufficiently regular in time, and the initial data are also sufficiently smooth, so time-singularities and singularities due to the smoothness of the initial data or the change of boundary conditions are not handled by the method.

References