# Broadband spectral analysis of Mrk 926 using multi-epoch X-ray observations

S. Chalise,1 A. M. Lohfink,1 J. Chauhan,1 T. D. Russell,2 D. J. K. Buisson3 L.
Mallick4
1Department of of Physics, Montana State University, P.O. Box 173840, Bozeman,
MT 59717-3840, USA
2INAF, Istituto di Astrofisica Spaziale e Fisica Cosmica, Via U. La Malfa 153,
I-90146 Palermo, Italy
3Department of Astronomy, University of Southampton, Southampton SO17 1BJ, UK
4Cahill Center for Astronomy and Astrophysics, California Institute of
Technology, Pasadena, CA 91125, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The X-ray spectra of some active galactic nuclei (AGN) show a soft X-ray
excess, emission in excess to the extrapolated primary X-ray continuum below 2
keV. Recent studies have shown that this soft excess can be described well as
originating from either a relativistic ionized reflection, the extreme
blurring of the reprocessed emission from the innermost region of the
accretion disk, or Comptonization from an optically thick and warm region
called the ‘warm corona’, in which electron scattering is the dominant source
of opacity. To constrain the origin of the soft excess in the Seyfert 1 galaxy
Mrk 926, we carry out an multi-epoch X-ray spectral study using observations
from Suzaku (2009), XMM-Newton and NuSTAR (2016), and NuSTAR and Swift-XRT
(2021). The broadband X-ray spectra of Mrk 926 contains: a thermally
Comptonized primary continuum, a variable soft excess, and distant reflection.
We find that in Mrk 926 as in so many sources, it is difficult to make a
definite statement as to what is causing the observed soft excess. A warm
coronal-like component is slightly preferred by the data but a reflection
origin is also possible. Using archival radio data, we detect an optically-
thin radio component in our broadband study of Mrk 926. While this component
is consistent with an optically-thin radio jet, future multi-wavelength
observations including high spatial resolution radio observations at multiple
frequencies are required to probe the origin of the radio emission in more
detail.

###### keywords:

Mrk 926 – MCG-02-28-22 – Warm corona – radiation mechanisms: non–thermal

††pubyear: 2022††pagerange: Broadband spectral analysis of Mrk 926 using
multi-epoch X-ray observations–References

##  1 Introduction

Active galactic nuclei (AGN) are powered by accretion onto the supermassive
black hole (SMBH) at the galactic center of their host galaxies. They are
highly energetic, and common X-ray sources in the universe. X-ray emission
from AGN, produced at the innermost region, is a great probe to study the
physical processes in the extreme gravity regime. A subgroup of AGN, the
Seyfert 1 galaxies, are particularly useful, in an observational sense, due to
the relatively unobscured view of their central engine.

All Seyfert 1 X-ray spectra share key similarities. Their primary X-ray
emission can be approximated as a cut-off powerlaw, and is consistent with
being produced via thermal Comptonization of optical-UV photons from the
accretion disk. This thermal Comptonization occurs in a compact region called
the hot corona (Haardt & Maraschi, 1991; Haardt et al., 1994). A part of this
primary emission can be reprocessed by interacting with different regions of
the accretion disk or the dusty torus, generating a ‘reflection’ spectrum.
Reflection from distant matter generates a neutral Fe K$\alpha$ fluorescent
emission line (George & Fabian, 1991; Matt et al., 1991) at 6.4 keV. If the
reflection is from an area closer to the SMBH, Doppler shift and gravitational
redshift broadens the fluorescent line profile (Fabian et al., 1989) and more
complex line emission is produced as the material is ionized. The reflection
spectrum is also accompanied by the “Compton hump" feature around 20–30 keV
due to Compton scattering of the primary continuum from the disk or parsec
(pc) scale torus (Ross et al., 1999).

In more than 50 $\%$ of Seyfert 1 galaxies, soft X-ray emission in excess to
that of the extrapolated hard X-ray continuum is observed below 2 keV and
referred to as ‘soft excess’. Its origin is still not completely understood,
and is currently theorized to be either due to warm Comptonization or
relativistic reflection from the accretion disk (Arnaud et al., 1985; Walter &
Fink, 1993; Piconcelli et al., 2005; Winter et al., 2012; Done et al., 2012).
The emitted X-rays can also be obscured by the gas and dust present in the
line-of-sight, imparting absorption features. All these components make up a
typical Seyfert X-ray spectrum. We can study the AGN phenomenon by analyzing
this spectrum using physically motivated models.

A group of nearby Seyfert 1 galaxies that show little to no X-ray absorption
are the so-called “bare Seyferts". Their study is free from the uncertainties
introduced from the complex modeling of the absorption providing a clear
measurement of the central region. Mrk 926, also known as MCG-2-58-22, is a
local (z=0.047) X-ray luminous “bare" AGN classified in the optical as a
Seyfert 1.5 galaxy. The lack of extended structures in the radio observations
of Mrk 926 (VLBA, 8 GHz) with a high temperature sub-pc scale core supports
its “bare" classification (Mundell et al., 2000). We usually expect cold and
diffuse radio emissions from dusty obscurers. Also, Mrk 926’s 2–10 keV X-ray
flux is known to vary on a timescale of months to years (Choi et al., 2001;
Choi et al., 2002). AGN variability can provide an additional probe into the
scale and structure of the innermost region. Thus, a clear view of a variable
central engine makes Mrk 926 an interesting X-ray source.

Mrk 926 has been observed previously with many X-ray observatories such as
ROSAT , ASCA , XMM , and BeppoSAX (Weaver et al., 1995; Choi et al., 2001;
Bianchi et al., 2004, etc.). Those observations have revealed the presence of
a variable soft excess (Ghosh & Soundararajaperumal, 1992), a variable Fe K
line profile (Weaver et al., 2001), and a Compton hump (Bianchi et al., 2004)
atop the primary continuum. This was also confirmed by the analysis of a 140
ks Suzaku observation (Rivers et al., 2011), which found a primary power-law
continuum, a soft excess and placed tight constraints on the strength and
nature of the reflection features. The authors also found no significant
evidence for broadening of the Fe line expected from the inner accretion disk
reflection (Rivers et al., 2011). The same Suzaku and the 2016 joint XMM
-NuSTAR observations were studied by Laha & Ghosh (2021), where the authors
found the vanishing of Compton-hump between 2009 and 2016 observations with
similar primary continuum. The authors proposed a scenario of a dynamic torus
structure as an explanation (Laha & Ghosh, 2021).

A recent velocity-resolved reverberation mapping of Mrk 926 was carried out by
Kollatschny et al. (2022) with a detailed spectroscopic and photometric
variability study of its very broad emission lines (FWHM $\gg 4000$ km s-1).
They found a drastic decrease of the optical continuum luminosity (50% in 2.5
months) suggesting high variability. The inclination angle of the line
emitting region was found to be $\sim$50∘. The black hole mass was derived to
be 1.1$(\pm{0.2})\times 10^{8}$M⊙, indicating a low Eddington ratio (3-8%).
Interestingly, they also found additional fast-response outer emission
components, the Balmer satellites, which originate in a different, spatially
distinct region such as a small-scale central radio jet.

There are still several open questions regarding Mrk 926’s X-ray spectrum,
mentioned in the following text. Using archival ROSAT and ASCA data, Choi et
al. (2001) found a broad iron line centered at $\sim$ 6.3 keV. This is not
unexpected from a bare nucleus as we have a clear view of the innermost
emission regions where the relativistic broadening of the iron line happens.
However, other studies such as Rivers et al. (2011) (140 ks Suzaku broadband
data) and Bianchi et al. (2004) (simultaneous $\sim 7$ ks XMM and BeppoSAX
broadband data) did not find a broadened line. Therefore, a pertinent question
is whether the broad iron line is a transient phenomenon in Mrk 926 or is its
strength variable and is sometimes overwhelmed by the rest of the emission? If
the iron line is variable, what is causing these changes? There is historical
evidence of strong soft X-ray variability in Mrk 926 (Ghosh &
Soundararajaperumal, 1992). The origin of the soft excess, and the physical
process responsible for this variability is still unknown. In this study of
Mrk 926, we make an effort to answer these open questions.

In this paper, we report the results from a 2009 Suzaku  observation, a 2016
overlapping XMM and NuSTAR campaign, and a 2021 NuSTAR and Swift -XRT
observation of Mrk 926. We analyze these multi-epoch spectra and shed more
light on the physical processes that govern the central engine of this AGN.
The next section tabulates the new observations and describes the data
reduction. We then perform a detailed spectral analysis in Section 3. We
subsequently discuss the results and their relevance for our understanding of
Mrk 926 and other Seyfert-like AGN in Section 4. Finally, we summarize our
findings in Section 5.

##  2 Observations and data reduction

###  2.1 Overview

In this paper, we analyze the broadband X-ray spectra of Mrk 926 observed in
three different epochs (Rivers et al., 2011; Laha & Ghosh, 2021, this work).
The first epoch is the 2009 Suzaku observation, the second epoch is the joint
2016 XMM -NuSTAR observation, and the third epoch is the 2021 NuSTAR (Harrison
et al., 2013) observation with a Swift -XRT pointing. The details of these
epochs are provided in table 1.

Table 1: The details of Mrk 926 X-ray observations used in this work. Satellite | ObsID | Date | Net Exposure | Epoch
---|---|---|---|---
Suzaku | 704032010 | 2009-12-02 | 139 ks | 1
XMM-Newton | 0790640101 | 2016-11-21 | 59 ks | 2
NuSTAR | 60201042002 | 2016-11-21 | 106 ks | 2
60761009002 | 2021-07-04 | 17 ks | 3
Swift-XRT | 00089294001 | 2021-07-04 | 2 ks | 3

###  2.2 Suzaku

The observation was performed by the XIS 0/1/3 CCD cameras and the Hard X-ray
Detector (HXD) Positive Intrinsic Negative (PIN) silicone diodes. Our data
reduction method is similar to Rivers et al. (2011), which also studied this
observation. We generated the cleaned event files with the Suzaku calibration
database released on 2018 October 10. The source spectrum was selected from a
circular region of 170 arcsec radius, while the background was from a circular
source-free region of 125 arcsec. The ARFs and RMFs were generated using
xisarfgen and xisrmfgen tasks, respectively. The HXD-PIN data were reduced
with aepipeline task. We utilized the "tuned" background event files to
produce the non X-ray background spectra of HXD-PIN, to which the simulated
cosmic X-ray background spectrum was added.

###  2.3 XMM-Newton

The XMM-Newton Observation data files (ODF) were processed using the XMM-
Newton Science Analysis System (XMM-SAS v20.0). The MOS detector, operated in
small window mode, is fully covered by the source photons and the selection of
a source-free background region is not feasible. We thus focus our analysis on
the EPIC-PN. The EPIC-PN data were first screened and periods of high particle
backgrounds rejected. The spectra and light curves were created from the
cleaned event files using evselect. We select an annular source region of 40
arcsec outer radius and 8 arcsec inner radius in order to mitigate mild pile-
up. We also select a circular background region of 60 arcsec. Responses were
created with arfgen and rmfgen tools. In arfgen, we apply the empirical
correction of the EPIC effective area with the parameter applyabsfluxcorr=yes.
This correction is essential to rectify the spectral shape mismatch seen with
the simultaneous NuSTAR observation.

###  2.4 NuSTAR

The NuSTAR data from both detectors, FPMA and FPMB, were reduced using the
standard pipeline ( NUPIPELINE ) of the NuSTAR Data Analysis Software (
NUSTARDAS v2.1.1). Calibration files from NuSTAR  CALDB  v20210202 were used.
We utilized the background filtering reports provided by the NuSTAR team to
reduce the effect of high particle background and produce cleaned event files.
We selected SAAMODE “optimized" and exclude the “tentacle" region minimizing
the effect of the South Atlantic Anomaly (SAA) and optimizing exposure time
and background level. We then extracted the spectra and light curves from a
circular source region of 40 arcsec and a circular background region of 100
arcsec close to the source.

###  2.5 Swift-XRT

The Swift -XRT spectra were obtained using the Swift -XRT data products
generator (Evans et al., 2009) supplied by the UK Swift Science Data Centre at
the University of Leicester.

##  3 Spectral analysis

In section 2, we generated the multi-epoch X-ray spectra of Mrk 926 using
Suzaku  (2009), joint XMM -NuSTAR  (2016) and joint NuSTAR -Swift  (2021)
observations. We bin them such that the detector resolution is over-sampled by
a factor of 3, ensuring the minimum signal to noise ratio of 10, except for
the Swift -XRT spectra which were binned to a minimum of 20 counts per bin. We
use the energy range of 0.85–1.7 keV, 2.3–10.0 keV for XIS; 13.0–45.0 keV for
PIN; 0.3–10.0 keV for PN; 3.0–75.0 keV for FPM; 0.3–10.0 keV for XRT. The
resulting spectra are plotted in the top subplot of Figure 1 with further
binning for clarity. We carry out the spectral analysis and model fitting with
XSPEC 12.11.1 (Arnaud & A., 1996) using the chi-square minimization technique;
all uncertainties are reported at 90% confidence.

In addition, our models always include: a cross correlation constant that
accounts for any cross-calibration flux offsets among the different spectra in
the same epoch; a galactic absorption component that accounts for the ISM
absorption from our galaxy using the ISM absorption model TBabs (Wilms et al.,
2000) with a column density ($N_{\text{H}}=2.86\times 10^{20}$ cm-2) obtained
from the N${}_{\text{H}}$ Tool (Kalberla et al., 2005) using the latest HI4PI
data. The cross-calibration constant for Suzaku -PIN is fixed at 1.16 which is
the expected value for XIS-nominal pointing.

A power-law spectrum is a rough description of an AGN primary continuum in the
X-ray band. Therefore to explore the additional X-ray components, We initially
fit a power-law to each epoch in the energy range 3–5 keV, where the spectra
are likely dominated by the primary continuum. We then apply this model,
without fitting, to the full energy range and search for residual features
(Figure 1, bottom). The three epochs show variations in both low and hard end
of the X-ray spectrum.

Figure 1: Top: All counts spectra for Suzaku, XMM, NuSTAR and Swift-XRT are
shown; the spectra have been binned further for plotting. The number in
brackets denote their respective epoch. Bottom: The ratio residuals after
fitting the spectra with a power-law in the range of 3.0–5.0 keV. The dashed
line indicates the neutral Fe-K$\alpha$ energy.

We notice smoothly-rising positive residuals towards the soft energies ($<$2
keV) indicative of a soft excess for Epoch 1 and 2. This could also be present
in Epoch 3 but simply not resolved by Swift-XRT. We also observe a line-like
feature at the neutral Fe-K$\alpha$ energy. At higher energies ($>$15 keV),
Epoch 2 shows a clear excess. These features can be associated with an X-ray
reflection spectrum, which are the reprocessed hard X-ray photons from the hot
corona after being backscattered off the AGN accretion disk or a distant
reflector. This reflection spectrum contains two main features: the 6.4 keV
Fe-K$\alpha$ line along with a bump-like feature that usually peaks around 25
keV. If this reflection occurs close to the black hole, these features can be
relativistically smeared. In some cases, this relativistic reflection can
explain the observed soft-excess (García et al., 2019). Thus, we check if the
reflection scenario can adequately explain the broadband multi-epoch spectra
of Mrk 926.

###  3.1 Dual reflection

AGNs are observed to be intrinsically variable in short timescales of a week
or less. So in each epoch, we let the primary continuum to be different.
However, we assume a stable reflecting medium between epochs as their
variability timescale is generally much larger.

We model the primary continuum using the thermal Comptonization model Nthcomp
(Życki et al., 1999). The powerlaw photon index ($\Gamma$) and its
normalization is allowed to vary among different epochs. The disk-blackbody
seed photon temperature (kTbb) is not sensitive to the fit and is fixed at
0.05 keV. The electron temperature (kTe) is a free parameter linked among
epochs. For reflection off this primary continuum, we use xillverCp  (Garcia
et al., 2013) for the distant reflection from the outer disk and torus, and
relxillCp  (Garcia et al., 2014) for relativistic reflection. The reflection
fraction for both are negative (reflection spectrum only) and variable among
epochs while the Power-law photon index ($\Gamma$), normalization and Electron
temperature (kTe) are linked to the primary continuum of the respective epoch.
The inclination angle (incl) of both reflectors is linked and only allowed to
vary between the range of 45∘–55∘. This inclination angle range was determined
by Kollatschny et al. (2022) using high-cadence spectroscopic variability
campaign of Mrk 926 in the optical band. The ionization parameter (logxi) of
the distant reflector is fixed at 0 while it is a free parameter for the
relativistic reflector. The iron abundance (A${}_{\text{fe}}$) for both
reflectors are free but not allowed to vary among epochs. For the relativistic
reflection: the inner (R${}_{\text{in}}$) and outer (R${}_{\text{out}}$)
radius of the accretion disk are fixed at the innermost stable circular orbit
(ISCO) and 400 R${}_{\text{g}}$ (gravitational radii) respectively; the
emissivity index up to 15 gravitational radii is a variable parameter among
epochs while the emissivity index beyond that is a free parameter linked among
epochs; the spin of the black hole (a) is fixed at the maximum value of 0.998.

Due to the limited photon counts, we find some reflector parameters to be
unconstrained for Epoch 3. Thus, we link all relxillCp  parameters except the
reflection fraction from Epoch 3 to Epoch 1. The fit for this dual reflection
scenario is satisfactory with a chi-square of 1582 for 1143 degrees of
freedom. We constrained the electron temperature (kTe) to $>284$ keV; the iron
abundance of the distant reflector to 5.0$\begin{subarray}{c}+2.3\\\
-0.3\end{subarray}$ times the solar value; the iron abundance of the
relativistic reflector to 1.6$\pm{0.4}$ times the solar value; the ionization
parameter (logxi) to 2.70$\pm{0.02}$. The detailed parameter values of the fit
are shown in Table 2. We note that the ratio residuals plot of this fit
(Figure 2, top) shows that it cannot explain the data above 30 keV perfectly.

###  3.2 Warm corona

The soft-excess can also be modeled as a Comptonized emission from a separate
component, often thought to be a warm and optically thick region called a warm
corona (Petrucci et al., 2018). We investigate if such a warm corona-like
component instead of a relativistic reflection can fit the multi-epoch
spectra.

We retain the same model structure for the primary continuum and its distant
reflection, and replace the relativistic reflection with a warm corona. We
model the warm corona using the thermal Comptonization model Nthcomp  (Życki
et al., 1999). The Power-law photon index ($\Gamma$), the electron temperature
(kTe) and the normalization for warm corona is allowed to vary among different
epochs. The disk-blackbody seed photon temperature (kTbb) is fixed at 0.05
keV. Similar to Section 3.1, we find the warm coronal parameters to be
unconstrained for Epoch 3. Thus, we link all except the normalization
parameter from Epoch 3 to Epoch 1. This model fits data well with a chi-square
of 1479 for 1139 degrees of freedom. We constrained the electron temperature
of the hot corona to 39$\pm{9}$ keV; the iron abundance of the non-
relativistic reflector to $>8.4$ times the solar value. The detailed parameter
values of the fit are shown in Table 3. The ratio residuals plot of this fit
(Figure 2, bottom) shows that this model explains the data above 30 keV better
than the dual reflection scenario. We note that this fit is slightly better
with the inclusion of additional narrow Gaussian lines in the XMM  spectra,
but for simplicity these were omitted as the results are unaffected.

Table 2: Spectral parameters for the dual reflection model (Section 3.1) containing the primary continuum, distant reflection and relativistic reflection fitted to the multi-epoch data. L1 denotes that the parameter is linked to the respective Epoch 1 parameter whereas $f$ means the parameter is fixed at the given value. Model | Parameter description | Symbol (Units) | Epoch 1 | Epoch 2 | Epoch 3
---|---|---|---|---|---
Primary Continuum (nthcomp) | Photon Index | $\Gamma$ | 1.83$\pm{0.06}$ | 1.88$\pm{0.05}$ | 1.73$\pm{0.04}$
Electron Temperature |  kTe (keV) |  $>$284 | L1 | L1
Normalization |  norm (10-3) | 14.1$\pm{0.2}$ | 12.2$\pm{0.2}$ | 4.8$\pm{0.3}$
Distant Reflection (xillverCp) | Iron Abundance |  Afe (solar) | 5.0$\begin{subarray}{c}+2.3\\\ -0.3\end{subarray}$ | L1 | L1
Ionization parameter |  logxi (ergs cm${}^{-}2$ s-1) | 0f | L1 | L1
Inclination |  incl (∘) | 53$\pm{2}$ | L1 | L1
Reflection fraction |  refl_frac (-10-3) | 3.0$\pm{0.4}$ | 5.5$\pm{0.6}$ | 7.4$\pm{3.5}$
Relativistic Reflection (relxillCp) | Inner emissivity index | Index1 | 4.5$\pm{0.3}$ |  $>$9.8 | L1
Outer emissivity index | Index2 |  $>$4.0 | L1 | L1
Inner-outer break |  Rbr (gravitational radii) | 15f | L1 | L1
Spin of the black hole | a | 0.998f | L1 | L1
Inner radius |  Rin (ISCO) | 1f | L1 | L1
Outer radius |  Rout (gravitational radii) | 400f | L1 | L1
Ionization parameter |  logxi (erg cm-2 s-1) | 2.70$\pm{0.02}$ | L1 | L1
Iron Abundance |  Afe (solar) | 1.6$\pm{0.4}$ | L1 | L1
Reflection fraction |  refl_frac (-10-3) | 6.1$\pm{0.5}$ | 8.4$\pm{0.3}$ |  $>$6.4
Figure 2: The top panel shows the ratio residuals for the dual reflection model fit in Section 3.1; the bottom panel shows the ratio residuals for the warm corona model fit in Section 3.2. The colors represent the same spectra as Figure 1.  Table 3: Spectral parameters for the warm corona model (Section 3.2) the primary continuum, distant reflection and a warm corona fitted to the multi-epoch data. L1 denotes that the parameter is linked to the respective Epoch 1 parameter whereas $f$ means the parameter is fixed at the given value. Model | Parameter description | Symbol (Units) | Epoch 1 | Epoch 2 | Epoch 3
---|---|---|---|---|---
Primary Continuum (nthcomp) | Photon Index | $\Gamma$ | 1.69$\pm{0.01}$ | 1.77$\pm{0.01}$ | 1.70$\pm{0.02}$
Electron Temperature |  kTe (keV) | 39$\pm{9}$ | L1 | L1
Normalization |  norm (10-3) | 12.7$\pm{0.5}$ | 13.1$\pm{0.1}$ | 4.8$\pm{0.2}$
Distant Reflection (xillverCp) | Iron Abundance |  Afe (solar) |  $>$8.4 | L1 | L1
Ionization parameter |  logxi (ergs cm${}^{-}2$ s-1) | 0f | L1 | L1
Inclination |  incl (∘) |  $<$55 | L1 | L1
Reflection fraction |  refl_frac (-10-3) | 2.6$\pm{0.4}$ | 2.3$\pm{0.4}$ | 4.6$\pm{2.2}$
Warm Corona (nthcomp) | Photon Index | $\Gamma$ | 2.84$\pm{0.1}$ | 2.63$\pm{0.1}$ | L1
Electron Temperature |  kTe (keV) |  $>$0.85 | 0.23$\pm{0.02}$ |  L1
Normalization |  norm (10-3) | 4.3$\pm{0.6}$ | 1.5$\pm{0.2}$ |  $<$0.16

###  3.3 Mrk 926 at Lower Frequencies

Using the Karl G. Jansky Very Large Array (VLA) observations at 4.8 GHz with
arcsec resolution, Ulvestad & Wilson (1984) found Mrk 926 to be a relatively
powerful radio source for a type 1 Seyfert with an unresolved core. Mundell et
al. (2000) studied the high angular resolution ($\sim$2 mas) radio continuum
observations of Mrk 926 using the Very Long Baseline Array (VLBA) at 8.4 GHz.
They found the core to be unresolved with brightness temperature in excess of
108 K and size less than 1 pc. They found the radio core to be consistent with
a non-thermal synchrotron self-absorption region without a Doppler boost. Mrk
926 was also monitored at 22 GHz with 1 arcsec resolution by the VLA (Smith et
al., 2020). Using these past radio observations at arcsec or sub-arcsec
resolution and additional archival flux measurements at other bands from Far-
IR up to X-rays, we study the nature of the self-absorbed synchrotron core.

We use the online Spectral SED builder tool111The SED builder is an online
tool dedicated to multi-frequency data visualization, together with fitting
routines useful for extracting refined scientific products. Provided by the
Space Science Data Center (SSDC): http://www.ssdc.asi.it to generate an
spectral energy distribution (SED) using archival flux measurements of Mrk
926. The SED is shown in Figure 3. We find the radio SED to be consistent with
a power-law with a steep spectral slope (F${}_{\nu}\propto\nu^{-0.54}$) which
is indicative of an optically-thin component (Eckart et al., 1986).

Figure 3: The SED generated by the online SED builder tool1 for Mrk 926 using
archival flux measurements. The blue line is a power-law fit to the arcsec or
higher resolution radio data.

In the optical reverberation mapping of Mrk 926 (Kollatschny et al., 2022), a
Balmer satellite component, consistent with the presence of a compact jet, was
found. A compact jet is continually-accelerated along the jet, where we detect
single synchrotron spectra from different distances along the jet resulting in
a flat to slightly inverted radio spectrum at GHz frequencies (Blandford &
Königl, 1979; Panessa et al., 2019). The flat/inverted radio spectrum will
break at some higher frequency $\nu_{b}$, either at the particle acceleration
region in a shocked zone or the base of the jet, into a steep spectrum
(Blandford & Königl, 1979; Koljonen et al., 2015).

There could be an underlying compact jet in Mrk 926 with the break frequency
present at MHz or higher frequencies. This jet structure, if present, could be
identified at higher spatial resolutions. However, with the available archival
data, we cannot postulate the origin of this radio emission. A future multi-
wavelength observing campaign of Mrk 926 including higher spatial resolution
radio observations at multiple frequencies is highly desired to study it in
detail. This will allow us to perform a full SED modeling including the disk,
host-galaxy etc. without the complications from variability, and resolve this
structure.

##  4 Discussion

###  4.1 Comparison with results from previous X-ray studies

Historically, Mrk 926 has been observed to show flux variations on longer
timescales (months and years). Using data from various X-ray satellites
covering more than 20 years (1977–1999), Choi et al. (2002) reported the X-ray
continuum flux variations in Mrk 926. We update their observed 2-10 keV X-ray
flux variations plot (see Figure 4). We find that Mrk 926 still shows
significant continuum variability.

Figure 4: Historic 2–10 keV flux variation plot adapted from Choi et al.
(2002) and updated with recent values.

Additionally, Mrk 926 also shows variations in its iron line profile (Weaver
et al., 2001; Bianchi et al., 2004). Choi et al. (2001), using ROSAT and ASCA
archival data, found the line profile dominated by a broad component. Bianchi
et al. (2004), using a $\sim$7 ks XMM and a simultaneous BeppoSAX observation,
found the presence of a Fe K$\alpha$ and a Fe XXVI line. Meanwhile Rivers et
al. (2011), using Suzaku observation, found the need for Fe K$\alpha$,
K$\beta$, and an Fe K shell absorption edge to explain the observed line
profile, and concluded that a broad component, if present, must be weak.

The 2009 Suzaku observation (Epoch 1) were described by Rivers et al. (2011)
using a torus model and an additional power-law component for soft excess. The
torus was found to have subsolar iron abundance with a column density,
log(N${}_{\texttt{H}}$) of 24.55$\pm{0.15}$ while the inclination angle and
the covering factor was fixed to 30∘ and 0.5 respectively. The same
observation was used by Noda et al. (2013) to study the nature of soft excess
in Mrk 926 using intensity-correlated spectral analysis. They found that the
thermal Comptonization interpretation can explain the broadband Suzaku data
including the soft excess, while the relativistic reflection interpretation
was failed to explain it. This is in agreement with our results in Section 3.1
and 3.2.

The 2009 Suzaku and 2016 XMM -NuSTAR  observations (Epoch 1 and 2) were
studied by Laha & Ghosh (2021). They modeled the soft excess with just the
relativistic reflection and found that between two epochs: the primary
continuum is constant both in photon index and luminosity; the soft excess
flux decreases by a factor of 2; the Compton hump vanishes; the narrow
Fe-k$\alpha$ line becomes marginally broad and its flux doubles. They argue
since the primary continuum does not change but the Compton hump vanishes
between two epochs, the reprocessing medium must be changing. They speculate
that this could be explained by a dynamic torus where the equatorial toroidal
disk produces the narrow Fe-k$\alpha$ emission while the outflowing component
generates the variable Compton hump.

To explore this hypothesis, we replace xillverCp  with a torus reprocessor
model borus12 (Baloković et al., 2018) in our multi-epoch analysis. Assuming
either a stable or a variable torus structure between epochs, we find minimal
changes to the fit and still see a small $>30$ keV excess in the dual
reflection scenario. For the warm corona scenario, we can get a fit similar to
our best fit with a stable torus reprocessor. Additionally, recent simulations
show that the outflowing component of the torus can change the equivalent
width of some X-ray fluorescent lines but does not significantly influence the
X-ray spectra (McKaig et al., 2022). Thus, we deem the scenario presented in
Laha & Ghosh (2021) to be unlikely considering the warm corona/dual reflection
model with a distant reflector, either torus or an accretion disk can describe
the multi-epoch broadband X-ray spectra of Mrk 926 without requiring
significant changes to the reprocessor between epochs.

###  4.2 Soft excess Origin

Recent studies debate two possible origins of AGN soft excess: the
relativistic ionized reflection of the continuum from the inner accretion
disk, and the Comptonization of the accretion disk photons in an optically
thick, warm corona (e.g. Petrucci et al., 2013; García et al., 2019).

There are many AGN for which the relativistic reflection origin of soft excess
works well (e.g. Mallick et al., 2018; Jiang et al., 2019; Middei et al.,
2020). In this scenario, the extreme gravitational blurring blends the rich
forest of fluorescent emission lines produced by the reprocessing of X-ray
photons in the innermost regions of accretion disk creating a smooth soft
excess. However, this scenario requires extreme parameters, such as maximum
black hole spin, a very low and compact hot corona, and a very high density
for the inner accretion disk (García et al., 2019). From Section 3, we find
that the relativistic reflection struggles to describe the broadband X-ray
spectra of Mrk 926 well above 30 keV. This could either mean that the soft
excess in Mrk 926 is not due to the relativistic reflection alone or that our
simplistic reflection modeling still has shortcomings.

Testing the warm coronal origin of soft excess in a sample of 22 unabsorbed
radio-quiet AGN with combined 100 observations, Petrucci et al. (2018) found
it to be a good description to more than 90% of the sample with the warm
coronal temperature in the range $\sim 0.1-1$ keV, and the optical depth in
$\sim 10-40$ range. However, at these temperatures and optical depth, the
atomic opacities is likely to dominate over the Thomson opacities and create
strong absorption features instead of a smooth soft excess (García et al.,
2019). This would also lead to rapid cooling of the warm corona, requiring a
well-tuned heating process that constantly offsets this (Ballantyne, 2020;
Ballantyne & Xiang, 2020). Petrucci et al. (2020) studied the radiative
equilibrium and emission of a warm corona in the presence of an internal
heating source and found that if this corona is covering a large part of a
weakly dissipative accretion disk, it can generate emission matching the
observed soft-excess. One such source of internal heating could be the
presence of magnetized accretion flow.

From Section 3, we find that the warm corona model can model the broadband
X-ray spectra of Mrk 926 in all three epochs well at all energies. Thus, it is
certainly a possibility for the soft excess origin in this source. To make a
strong statement about the nature of the soft excess in this source, high
quality data covering the 50-100 keV range would be needed. Finally, we note
that since we are unable to constrain the nature of radio emitter at high
frequencies, we cannot fully exclude any flux contribution at soft X-rays.

##  5 Summary

We present results of a broadband X-ray study of Mrk 926 using the 2009 Suzaku
(Epoch 1), 2016 XMM -NuSTAR  (Epoch 2), and 2021 NuSTAR -Swift -XRT (Epoch 3)
observations. Our findings can be summarized as follows:

  1. 1.

Mrk 926 is a variable AGN. The unabsorbed 2-10 keV luminosity of Mrk 926 is
2.9$\pm{0.1}\times 10^{44}$ erg s-1 (Epoch 1 and 2) and 1.1$\pm{0.1}\times
10^{44}$ erg s-1 (Epoch 3).

  2. 2.

We found the broadband spectra of Mrk 926 show a primary continuum, a soft
excess and distant reflection. If the soft excess is modeled using
relativistic reflection, we find some residuals above 30 keV. We find no such
residuals if we describe soft excess as a thermally Comptonized emission from
a optically-thick warm corona.

  3. 3.

An optical Balmer satellite component was found in previous study of Mrk 926
which suggests the presence of a compact jet. The signature of a compact jet
in radio band is the presence of a flat/inverted spectrum that breaks at some
higher frequency to a steep spectrum. We found a steep radio spectrum in the
archival SED of Mrk 926. Thus, the archival radio data cannot confirm the
presence of a compact jet.

Future multi-wavelength observing campaign including a high spatial resolution
radio observation is highly desired to study the origin of this radio emitter
in the ‘bare’ Seyfert Mrk 926.

## Acknowledgements

This work made use of data from the NuSTAR mission, a project led by the
California Institute of Technology, managed by the Jet Propulsion Laboratory,
and funded by the National Aeronautics and Space Administration. This research
has also made use of the NASA/IPAC Extragalactic Database (NED), which is
funded by the National Aeronautics and Space Administration and operated by
the California Institute of Technology. We would like to thank Javier García
and Andy Fabian for helpful discussions.

## Data Availability

All the data used in this article are publicly available from ESA XMM-Newton
Science Archive (XSA; http://nxsa.esac.esa.int/) and NASA High Energy
Astrophysics Science Archive Research Center (HEASARC;
https://heasarc.gsfc.nasa.gov/). The reduced data products used in this work
may be shared on reasonable request to the authors.

## References

  * Arnaud & A. (1996) Arnaud A. K., 1996, Astronomical Data Analysis Software and Systems V, 101, 17
  * Arnaud et al. (1985) Arnaud K. A., et al., 1985, MNRAS, 217, 105
  * Ballantyne (2020) Ballantyne D. R., 2020, MNRAS, 491, 3553
  * Ballantyne & Xiang (2020) Ballantyne D. R., Xiang X., 2020, MNRAS, 496, 4255
  * Baloković et al. (2018) Baloković M., et al., 2018, ApJ, 854, 42
  * Bianchi et al. (2004) Bianchi S., Matt G., Balestra I., Guainazzi M., Perola G. C., 2004, A&A, 422, 65
  * Blandford & Königl (1979) Blandford R. D., Königl A., 1979, ApJ, 232, 34
  * Choi et al. (2001) Choi C.-S., Dotani T., Yi I., Fletcher A., Kim C., 2001, Journal of Korean Astronomical Society, 34, 129
  * Choi et al. (2002) Choi C.-S., Dotani T., Chang H.-Y., Yi I., 2002, Journal of Korean Astronomical Society, 35, 1
  * Done et al. (2012) Done C., Davis S. W., Jin C., Blaes O., Ward M., 2012, MNRAS, 420, 1848
  * Eckart et al. (1986) Eckart A., Witzel A., Biermann P., Johnston K. J., Simon R., Schalinski C., Kuhr H., 1986, A&A, 168, 17
  * Evans et al. (2009) Evans P. A., et al., 2009, MNRAS, 397, 1177
  * Fabian et al. (1989) Fabian A. C., Rees M. J., Stella L., White N. E., 1989, MNRAS, 238, 729
  * Garcia et al. (2013) Garcia J., Dauser T., Reynolds C. S., Kallman T. R., McClintock J. E., Wilms J., Eikmann W., 2013, The Astrophysical Journal, 768, 146
  * Garcia et al. (2014) Garcia J., et al., 2014, The Astrophysical Journal, 782, 76
  * García et al. (2019) García J. A., et al., 2019, ApJ, 871, 88
  * George & Fabian (1991) George I. M., Fabian A. C., 1991, MNRAS, 249, 352
  * Ghosh & Soundararajaperumal (1992) Ghosh K. K., Soundararajaperumal S., 1992, ApJ, 398, 501
  * Haardt & Maraschi (1991) Haardt F., Maraschi L., 1991, ApJ, 380, L51
  * Haardt et al. (1994) Haardt F., Maraschi L., Ghisellini G., 1994, ApJ, 432, L95
  * Harrison et al. (2013) Harrison F. A., et al., 2013, The Astrophysical Journal, 770, 103
  * Jiang et al. (2019) Jiang J., et al., 2019, MNRAS, 489, 3436
  * Kalberla et al. (2005) Kalberla P. M. W., et al., 2005, Astronomy & Astrophysics, 440, 775
  * Koljonen et al. (2015) Koljonen K. I. I., et al., 2015, ApJ, 814, 139
  * Kollatschny et al. (2022) Kollatschny W., et al., 2022, A&A, 657, A122
  * Laha & Ghosh (2021) Laha S., Ghosh R., 2021, ApJ, 915, 93
  * Mallick et al. (2018) Mallick L., et al., 2018, MNRAS, 479, 615
  * Matt et al. (1991) Matt G., Perola G. C., Piro L., 1991, A&A, 247, 25
  * McKaig et al. (2022) McKaig J., Ricci C., Paltani S., Satyapal S., 2022, MNRAS, 512, 2961
  * Middei et al. (2020) Middei R., et al., 2020, A&A, 640, A99
  * Mundell et al. (2000) Mundell C. G., Wilson A. S., Ulvestad J. S., Roy A. L., 2000, ApJ, 529, 816
  * Noda et al. (2013) Noda H., Makishima K., Nakazawa K., Uchiyama H., Yamada S., Sakurai S., 2013, PASJ, 65, 4
  * Panessa et al. (2019) Panessa F., Baldi R. D., Laor A., Padovani P., Behar E., McHardy I., 2019, Nature Astronomy, 3, 387
  * Petrucci et al. (2013) Petrucci P. O., et al., 2013, A&A, 549, A73
  * Petrucci et al. (2018) Petrucci P. O., Ursini F., De Rosa A., Bianchi S., Cappi M., Matt G., Dadina M., Malzac J., 2018, A&A, 611, A59
  * Petrucci et al. (2020) Petrucci P. O., et al., 2020, A&A, 634, A85
  * Piconcelli et al. (2005) Piconcelli E., Jimenez-Bailón E., Guainazzi M., Schartel N., Rodríguez-Pascual P. M., Santos-Lleó M., 2005, A&A, 432, 15
  * Rivers et al. (2011) Rivers E., Markowitz A., Rothschild R., 2011, ApJ, 732, 36
  * Ross et al. (1999) Ross R. R., Fabian A. C., Young A. J., 1999, MNRAS, 306, 461
  * Smith et al. (2020) Smith K. L., et al., 2020, MNRAS, 492, 4216
  * Ulvestad & Wilson (1984) Ulvestad J. S., Wilson A. S., 1984, ApJ, 278, 544
  * Walter & Fink (1993) Walter R., Fink H. H., 1993, A&A, 274, 105
  * Weaver et al. (1995) Weaver K. A., Nousek J., Yaqoob T., Hayashida K., Murakami S., 1995, ApJ, 451, 147
  * Weaver et al. (2001) Weaver K. A., Gelbord J., Yaqoob T., 2001, ApJ, 550, 261
  * Wilms et al. (2000) Wilms J., Allen A., McCray R., 2000, The Astrophysical Journal, 542, 914
  * Winter et al. (2012) Winter L. M., Veilleux S., McKernan B., Kallman T. R., 2012, ApJ, 745, 107
  * Życki et al. (1999) Życki P. T., Done C., Smith D. A., 1999, MNRAS, 309, 561


# A SDFEM for system of two singularly perturbed problems of convection-
diffusion type with discontinuous source term.

A. Ramesh Babu
Department of Mathematics
School of Engineering
AMRITA Vishwa Vidyapeetham
Coimbatore - 641 112
Tamilnadu, India

###### Abstract.

We consider a system of two singularly perturbed Boundary Value Problems
(BVPs) of convection-diffusion type with discontinuous source terms and a
small positive parameter multiplying the highest derivatives. Then their
solutions exhibit boundary layers as well as weak interior layers. A numerical
method based on finite element method (Shishkin and Bakhvalov-Shishkin meshes)
is presented. We derive an error estimate of order $O(N^{-1}\ln^{3/2}{N})$ in
the energy norm with respect to the perturbation parameter. Numerical
experiments are also presented to support our theoritical results.
AMS Mathematics Subject Classification: 65L10, CR G1.7

Key words:Singularly perturbed problem, Discontinuous source term, Weakly
coupled system, Finite element method, Energy norm, Convection-diffusion,
Boundary value problem.

The author wishes to acknowledge a support rendered by Bharathidasan
University, Tamilnadu, India through collabration with Prof. N. Ramanujam.
E-mail<EMAIL_ADDRESS>matramesh2k5@yahoo.co.in(A. Ramesh Babu)

##  1\. Introduction

Singularly Perturbed Differential Equations(SPDEs) appear in several branches
of applied mathematics. Analytical and numerical treatment of these equations
have drawn much attention of many researchers [1, 3, 2, 4, 5]. In general,
classical numerical methods fail to produce good approximations for these
equations. Hence one has to look for non-classical methods. A good number of
articles have been appearing in the past three decades on non-classical
methods which cover mostly second order equations. But only a few authors have
developed numerical methods for singularly perturbed system of ordinary
differential equations.[7, 8, 10, 11, 12, 13].
Systems of this kind have applications in electro analytic chemistry when
investigating diffusion processes complicated by chemical reactions. The
parameters multiplying the highest derivatives characterize the diffusion
coefficient of the substances. Other applications include equations of
predator-prey population dynamics. As was mentioned above, classical numerical
methods fails to produce good approximations for singularly perturbed system
of equations also. Hence various methods are proposed in the literature in
order to obtain numerical solution to singularly perturbed system of second
order differential equations subject to Dirichlet type boundary conditions
when the source terms are smooth on $(0,1)$ [8, 11, 12]. Motivated by the
works of T. Lin$\ss{}$ and N. Madden [7], in the present paper we suggest a
numerical method for singularly perturbed weakly coupled system of two
ordinary differential equations of convection-diffusion type with
discontinuous source terms. Basically the method is based on Streamline
Diffusion Finite Element Method (SDFEM) with layer adapted meshes like
Shishkin and Bakhvalov-Shishkin meshes. For this method we derive an error
estimate of order $O(N^{-1}\ln^{3/2}{N})$ in the energy norm.
In this paper, we consider the system of singularly perturbed BVP with
discontinuous source term

$\displaystyle P_{1}\bar{u}:=-\varepsilon
u_{1}^{\prime\prime}(x)+b_{1}(x)u_{1}^{\prime}(x)+a_{11}(x)u_{1}(x)+a_{12}(x)u_{2}(x)=f_{1}(x),\quad
x\in(\Omega^{-}\cup\Omega^{+})$ (1.1) $\displaystyle
P_{2}\bar{u}:=-\varepsilon
u_{2}^{\prime\prime}(x)+b_{2}(x)u_{2}^{\prime}(x)+a_{21}(x)u_{1}(x)+a_{22}(x)u_{2}(x)=f_{2}(x),\quad
x\in(\Omega^{-}\cup\Omega^{+})$ (1.2) $\displaystyle u_{1}(0)=0,\quad
u_{1}(1)=0,\quad u_{2}(0)=0,\quad u_{2}(1)=0,$ (1.3)

with the following conditions.

$\displaystyle b_{1}(x)\geq\beta_{1}>0,\quad b_{2}(x)\geq\beta_{2}>0,$ (1.4)
$\displaystyle a_{12}(x)\leq 0,\quad a_{21}(x)\leq 0,$ (1.5) $\displaystyle
a_{11}(x)>|a_{21}(x)|,\quad a_{22}(x)>|a_{12}(x)|,\quad\forall
x\in\bar{\Omega},$ (1.6)

$A=[a_{ij}],i=1,2;j=1,2\,$ satisfies the property

$\xi^{T}A\xi\geq\alpha\xi\xi^{T}\quad\text{for
every}\quad\xi=(\xi_{1},\xi_{2})\in\Re^{2}.$ (1.7)

For $k=1,2$

$\alpha-\frac{1}{2}b_{k}^{\prime}\geq\sigma_{k},\quad\text{for
some}\quad\alpha,\sigma_{k}>0.$ (1.8)

where $\varepsilon>0$ is a small parameter, $\Omega=(0,1),$
$\Omega^{-}=(0,d),$ $\Omega^{+}=(d,1),$ $d\in\Omega,$ and $u_{1},u_{2}\in
U\equiv C^{0}(\bar{\Omega})\cap C^{1}(\Omega)\cap
C^{2}(\Omega^{-}\cup\Omega^{+}),$ $\bar{u}=(u_{1},u_{2})^{T}$. Further it is
assumed that the source terms $f_{1},f_{2}$ are sufficiently smooth on
$\bar{\Omega}\setminus\\{d\\};$ both the functions $f_{1}(x)$ and $f_{2}(x)$
are assumed to have a single discontinuity at the point $d\in\Omega.$ In
general this discontinuity gives rise to interior layers in the solution of
the problem. Because $f_{i},i=1,2$ are discontinuous at $d$ the solution
$\bar{u}$ of (1.1) - (1.3) does not necessarily have a continuous second
derivative at the point $d.$ That is $u_{1},u_{2}\notin C^{2}(\Omega).$ But
the first derivative of the solution exists and is continuous. The authors
from [13] proved almost first order of convergence with respect to
$\varepsilon$ on a Shishkin mesh of the finite difference method with special
discretization in the point $d.$

######  Remark 1.1.

Through out this paper, $C,\,C_{1}$ denote generic constants that are
independent of the parameter $\varepsilon$ and $N,$ the dimension of the
discrete problem. We also assume $\varepsilon\leq CN^{-1}$ as is generally the
case in practice.

For our later analysis it is useful to have a decomposition of $\bar{u}$ in
the smooth part $\bar{v}$ and the layer part $\bar{w}.$ That is

$\bar{u}=\bar{v}+\bar{w}_{1}+\bar{w}_{2},\quad\text{where}\quad\bar{v}=(v_{1},v_{2}),\quad\bar{w}_{1}=(w_{11},w_{12}),\quad\bar{w}_{2}=(w_{21},w_{22}).$

######  Theorem 1.2.

With the decomposition of the above, for each $k,$ $0\leq k\leq 3,$ and
$j=1,2$ it holds

$\displaystyle|v_{j}^{(k)}(x)|\leq C(1+\varepsilon^{(2-k)}),\quad x\in\Omega,$
$\displaystyle|w_{1j}^{(k)}(x)|\leq
C\varepsilon^{-k}e^{\frac{-\beta(1-x)}{\varepsilon}},\quad x\in\bar{\Omega},$
$\displaystyle|w_{2j}^{(k)}(x)|\leq\begin{cases}C\varepsilon^{(1-k)}e^{\frac{-\beta(d-x)}{\varepsilon}},\quad
x\in\Omega^{-},\\\
C\varepsilon^{(1-k)}e^{\frac{-\beta(1-x)}{\varepsilon}},\quad
x\in\Omega^{+},\end{cases}$

where $\beta=\min\\{\beta_{1},\beta_{2}\\}.$

###### Proof.

Using the results of [10] and adopting the technique of [1] this theorem can
be proved. ∎

This paper is organized as follows. Section $2$ presents a weak formulation of
the BVP (1.1) - (1.3). We define an energy norm on $(H_{0}^{1}(\Omega))^{2}$
and discribe a finite element discretization of the problem. Section $3$
presents an analysis of the corresponding scheme on Shishkin and Bakhvalov-
Shishkin meshes. In section $4,$ we present an interpolation error on various
norms. The paper concludes with numerical examples.

##  2\. Analytical results

A standard weak formulation of (1.1)-(1.3) is: Find $u_{1},u_{2}\in
H_{0}^{1}(\Omega)$ such that

$\displaystyle B_{1}(u_{1},v_{1})=f_{1}(v_{1}),\quad\forall v_{1}\in
H_{0}^{1}(\Omega)$ (2.1) $\displaystyle
B_{2}(u_{2},v_{2})=f_{2}(v_{2}),\quad\forall v_{2}\in H_{0}^{1}(\Omega)$ (2.2)

where

$\displaystyle
B_{1}(u_{1},v_{1}):=\varepsilon(u_{1}^{\prime},v_{1}^{\prime})+(b_{1}u_{1}^{\prime},v_{1})+(a_{11}u_{1}+a_{12}u_{2},v_{1}),$
$\displaystyle
B_{2}(u_{2},v_{2}):=\varepsilon(u_{2}^{\prime},v_{2}^{\prime})+(b_{2}u_{2}^{\prime},v_{2})+(a_{21}u_{1}+a_{22}u_{2},v_{2})$

and

$\displaystyle f_{1}(v_{1})=(f_{1},v_{1}),$ $\displaystyle
f_{2}(v_{2})=(f_{2},v_{2}).$

Here $H_{0}^{1}(\Omega)$ denotes the usual Sobolev space and $(.,.)$ is the
inner product on $L_{2}(\Omega).$ Now we combine the two equations (2.1) -
(2.2) and get a single weak formulation. Then our problem is: Find
$\bar{u}\in(H_{0}^{1}(\Omega))^{2}$ such that

$B(\bar{u},\bar{v})=f(\bar{v}),\quad\forall\bar{v}\in(H_{0}^{1}(\Omega))^{2}$
(2.3)

with $B(\bar{u},\bar{v}):=B_{1}(u_{1},v_{1})+B_{2}(u_{2},v_{2})\quad$ and
$\quad f(\bar{v}):=f_{1}(v_{1})+f_{2}(v_{2}).$ Now we define a norm on
$(H_{0}^{1}(\Omega))^{2}$ associated with the bilinear form $B(.,.)$, called
continuous energy norm as
$|||\bar{u}|||_{H_{0}^{1}}=[\varepsilon(|u_{1}|_{1}^{2}+|u_{2}|_{1}^{2})+\sigma(\|u_{1}\|_{0}^{2}+\|u_{2}\|_{0}^{2})]^{1/2},$
where $\sigma=\min\\{\sigma_{1},\sigma_{2}\\}$ and $\|u\|_{0}:=(u,u)^{1/2}$ is
the standard norm on $L_{2}(\Omega),$ while $|u|_{1}:=\|u^{\prime}\|_{0}$ is
the usual semi-norm on $H_{0}^{1}(\Omega).$ We also use the notation
$\|\bar{u}\|_{0}=(\|u_{1}\|_{0}^{2}+\|u_{2}\|_{0}^{2})^{1/2}$ for the norm in
$(L_{2}(\Omega))^{2}.$
$B$ is a bilinear functional defined on $(H_{0}^{1}(\Omega))^{2}.$ Further we
have to prove that it is coercive with respect to $|||.|||_{H_{0}^{1}},$ that
is $B(\bar{u},\bar{u})\geq|||\bar{u}|||_{H_{0}^{1}}^{2}.$

######  Lemma 2.1.

A bilinear functional $B$ satisfies the coercive property with respect to
$|||.|||_{H_{0}^{1}}.$

###### Proof.

Let $\bar{u}=(u_{1},u_{2})\in(H_{0}^{1}(\Omega))^{2}.$ Then

$\displaystyle B(\bar{u},\bar{u})$ $\displaystyle=$
$\displaystyle\varepsilon(u_{1}^{\prime},u_{1}^{\prime})+(b_{1}u_{1}^{\prime},u_{1})+(a_{11}u_{1}+a_{12}u_{2},u_{1})+\varepsilon(u_{2}^{\prime},u_{2}^{\prime})+(b_{2}u_{2}^{\prime},u_{2})$
$\displaystyle+(a_{21}u_{1}+a_{22}u_{2},u_{2})$ $\displaystyle\geq$
$\displaystyle\varepsilon(|u_{1}|_{1}^{2}+|u_{2}|_{1}^{2})+\int_{0}^{1}b_{1}(x)u_{1}^{\prime}u_{1}dx+\int_{0}^{1}b_{2}(x)u_{2}^{\prime}u_{2}dx+(\alpha
u_{1},u_{1})$ $\displaystyle+(\alpha u_{2},u_{2})$ $\displaystyle=$
$\displaystyle\varepsilon(|u_{1}|_{1}^{2}+|u_{2}|_{1}^{2})+\int_{0}^{1}\frac{b_{1}(x)}{2}\frac{d}{dx}(u_{1}^{2})+\int_{0}^{1}\alpha
u_{1}^{2}dx+\int_{0}^{1}\frac{b_{2}(x)}{2}\frac{d}{dx}(u_{2}^{2})$
$\displaystyle+\int_{0}^{1}\alpha u_{2}^{2}dx$ $\displaystyle=$
$\displaystyle\varepsilon(|u_{1}|_{1}^{2}+|u_{2}|_{1}^{2})-\frac{1}{2}\int_{0}^{1}u_{1}^{2}d(b_{1}(x))+\int_{0}^{1}\alpha
u_{1}^{2}dx-\frac{1}{2}\int_{0}^{1}u_{2}^{2}d(b_{2}(x))$
$\displaystyle+\int_{0}^{1}\alpha u_{2}^{2}dx$ $\displaystyle=$
$\displaystyle\varepsilon(|u_{1}|_{1}^{2}+|u_{2}|_{1}^{2})+\int_{0}^{1}(\alpha-\frac{1}{2}b_{1}^{\prime}(x))u_{1}^{2}dx+\int_{0}^{1}(\alpha-\frac{1}{2}b_{2}^{\prime}(x))u_{2}^{2}dx$
$\displaystyle\geq$
$\displaystyle\varepsilon(|u_{1}|_{1}^{2}+|u_{2}|_{1}^{2})+\min\\{\sigma_{1},\sigma_{2}\\}[\int_{0}^{1}u_{1}^{2}dx+\int_{0}^{1}u_{1}^{2}dx]$
$\displaystyle B(\bar{u},\bar{u})$ $\displaystyle\geq$
$\displaystyle\varepsilon(|u_{1}|_{1}^{2}+|u_{2}|_{1}^{2})+\sigma(\|u_{1}\|_{0}^{2}+\|u_{2}\|_{0}^{2})$

Therefore we have

$B(\bar{u},\bar{u})\geq|||\bar{u}|||^{2}.$

Hence $B$ is coercive with respect to $|||.|||.$ ∎

Also $B$ is continuous in the energy norm and $f$ is a bounded linear
functional on $(H_{0}^{1}(\Omega))^{2}.$ By Lax-Milgram Theorem, we conclude
that the problem (2.3) has a unique solution.

###  2.1. Discretization of weak problem

Let $\Omega_{\varepsilon}^{N}=\\{x_{0},x_{1},\cdots,x_{N}\\}$ to be the set of
mesh points $x_{i}$, for some positive integer $N$. For
$i\in\\{1,2,\cdots,N\\}.$ We set $h_{i}=x_{i}-x_{i-1}$ to be the local mesh
step size, and for $i\in\\{1,2,\cdots,N\\}$ let
$\bar{h_{i}}=(h_{i}+h_{i+1})/2$. Let $V_{h}\subset H_{0}^{1}(\Omega)$ be the
space of piecewise linear functions on $\Omega$. As usual, basis functions of
$V_{h}$ are given by

$\phi_{i}(x)=\begin{cases}\frac{x-x_{i-1}}{h_{i}},\quad x\in[x_{i-1},x_{i}]\\\
\frac{x_{i+1}-x}{h_{i+1}},\quad x\in[x_{i},x_{i+1}]\\\ 0,\quad
x\notin[x_{i-1},x_{i+1}].\end{cases}$

Then our discretization of (2.3) is: Find $\bar{u}_{h}\in V_{h}^{2}$ such that

$B_{h}(\bar{u}_{h},\bar{v}_{h})=f_{h}(\bar{v}_{h}),\quad\forall\bar{v}_{h}\in
V_{h}^{2},$ (2.4)

where

$\displaystyle B_{h}(\bar{u}_{h},\bar{v}_{h}):=$ $\displaystyle(\varepsilon
u_{1h}^{\prime},v_{1h}^{\prime})+(b_{1}u_{1h}^{\prime},v_{1h})+(a_{11}u_{1h}+a_{12}u_{2h},v_{1h})+(\varepsilon
u_{2h}^{\prime},v_{2h}^{\prime})$
$\displaystyle+(b_{2}u_{2h}^{\prime},v_{2h})+(a_{21}u_{1h}+a_{22}u_{2h},v_{2h})$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(-\varepsilon
u_{1h}^{\prime\prime}(x)+b_{1}(x)u_{1h}^{\prime}(x)+a_{11}(x)u_{1h}(x)+a_{12}(x)u_{2h}(x))b_{1}v_{1h}^{\prime}dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(-\varepsilon
u_{2h}^{\prime\prime}(x)+b_{2}(x)u_{2h}^{\prime}(x)+a_{21}(x)u_{1h}(x)+a_{22}(x)u_{2h}(x))b_{2}v_{2h}^{\prime}dx$
$\displaystyle f_{h}(\bar{v}_{h}):=$
$\displaystyle(f_{1},v_{1h})+(f_{2},v_{2h})+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}f_{1}b_{1}v_{1h}^{\prime}+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}f_{2}b_{2}v_{2h}^{\prime}dx.$

The parameters $\delta_{1,i}\geq 0$ and $\delta_{2,i}\geq 0$ are called the
streamline-diffusion parameters and will be determined later. Here we define a
discrete energy norm on $V_{h}^{2}$ associated with the bilinear form
$B_{h}(.,.)$ as

$\displaystyle|||\bar{u}_{h}|||_{V_{h}}$ $\displaystyle=$
$\displaystyle[\varepsilon(|u_{1h}|_{1}^{2}+|u_{2h}|_{1}^{2})+\sigma(\|u_{1h}\|_{0}^{2}+\|u_{2h}\|_{0}^{2})+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}b_{1}^{2}(x_{i})(u_{1h}^{\prime}(x))^{2}dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}b_{2}^{2}(x_{i})(u_{2h}^{\prime}(x))^{2}dx]^{1/2}.$

$B_{h}$ is a bilinear functional defined on $V_{h}^{2}.$ Further we have to
prove that it is coercive with respect to $|||.|||_{V_{h}},$ that is
$B_{h}(\bar{u}_{h},\bar{u}_{h})\geq|||\bar{u}_{h}|||_{V_{h}}^{2}.$

######  Lemma 2.2.

If $\quad\delta_{1,i}=\delta_{2,i}=0\,\,\quad$ then $\quad
B_{h}(\bar{u}_{h},\bar{u}_{h})\geq|||\bar{u}_{h}|||_{V_{h}}^{2}$
and if $\quad
0<\delta_{1,i},\delta_{2,i}\leq\frac{1}{4}\min_{i=1,2}\\{\frac{\sigma_{i}}{\mu^{2}}\\},\quad\mu=\max_{x\in\bar{\Omega}}\\{\mid
a_{ij}(x)\mid\\},i,j=1,2\quad$ then $\quad
B_{h}(\bar{u}_{h},\bar{u}_{h})\geq\frac{1}{2}|||\bar{u}_{h}|||_{V_{h}}^{2}.$
That is, a bilinear functional $B_{h}$ satisfies the coercive property with
respect to $|||.|||_{V_{h}}.$

###### Proof.

Let $\bar{u}_{h}=(u_{1h},u_{2h})\in V_{h}^{2}.$ If
$\quad\delta_{1,i}=\delta_{2,i}=0\quad$ then the result directly follows from
Lemma (2.1).
If $\quad
0<\delta_{1,i},\delta_{2,i}\leq\frac{1}{4}\min_{i=1,2}\\{\frac{\sigma_{i}}{\mu^{2}}\\}\quad$
then we have

$\displaystyle B_{h}(\bar{u}_{h},\bar{u}_{h})$ $\displaystyle=$
$\displaystyle\varepsilon(u_{1h}^{\prime},u_{1h}^{\prime})+(b_{1}u_{1h}^{\prime},u_{1h})+(a_{11}u_{1h}+a_{12}u_{2h},u_{1h})+\varepsilon(u_{2h}^{\prime},u_{2h}^{\prime})+(b_{2}u_{2h}^{\prime},u_{2h})$
$\displaystyle+(a_{21}u_{1h}+a_{22}u_{2h},u_{2h})+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(-\varepsilon
u_{1h}^{\prime\prime}+b_{1}u_{1h}^{\prime}+a_{11}u_{1h}+a_{12}u_{2h})b_{1}u_{1h}^{\prime}dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(-\varepsilon
u_{2h}^{\prime\prime}+b_{2}u_{2h}^{\prime}+a_{21}u_{1h}+a_{22}u_{2h})b_{2}u_{2h}^{\prime}dx$
$\displaystyle\geq$
$\displaystyle\varepsilon(|u_{1h}|_{1}^{2}+|u_{2h}|_{1}^{2})+\int_{0}^{1}b_{1}(x)u_{1h}^{\prime}u_{1h}dx+\int_{0}^{1}b_{2}(x)u_{2h}^{\prime}u_{2h}dx+\int_{0}^{1}\alpha
u_{1h}^{2}dx$ $\displaystyle+\int_{0}^{1}\alpha
u_{2h}^{2}dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(b_{1}^{2}(u_{1h}^{\prime})^{2})dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(a_{11}u_{1h}+a_{12}u_{2h})b_{1}u_{1h}^{\prime}dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(b_{2}^{2}(u_{2h}^{\prime})^{2})dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(a_{21}u_{1h}+a_{22}u_{2h})b_{2}u_{2h}^{\prime}dx$
$\displaystyle\geq$
$\displaystyle\varepsilon(|u_{1h}|_{1}^{2}+|u_{2h}|_{1}^{2})+\int_{0}^{1}(\alpha-\frac{1}{2}b_{1}^{\prime}(x))u_{1h}^{2}dx+\int_{0}^{1}(\alpha-\frac{1}{2}b_{2}^{\prime}(x))u_{2h}^{2}dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(b_{1}^{2}(u_{1h}^{\prime})^{2})dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(b_{2}^{2}(u_{2h}^{\prime})^{2})dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(a_{11}u_{1h}+a_{12}u_{2h})b_{1}u_{1h}^{\prime}dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(a_{21}u_{1h}+a_{22}u_{2h})b_{2}u_{2h}^{\prime}dx$
$\displaystyle B_{h}(\bar{u}_{h},\bar{u}_{h})$ $\displaystyle\geq$
$\displaystyle\varepsilon(|u_{1h}|_{1}^{2}+|u_{2h}|_{1}^{2})+\sigma(\|u_{1h}\|_{0}^{2}+\|u_{2h}\|_{0}^{2})$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(b_{1}^{2}(u_{1h}^{\prime})^{2})dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(b_{2}^{2}(u_{2h}^{\prime})^{2})dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(a_{11}u_{1h}+a_{12}u_{2h})b_{1}u_{1h}^{\prime}dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(a_{21}u_{1h}+a_{22}u_{2h})b_{2}u_{2h}^{\prime}dx$

Using the assumption on $\delta_{1,i}$ and $\delta_{2,i},$ we obtain

$\displaystyle\mid\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(a_{11}u_{1h}+a_{12}u_{2h})b_{1}u_{1h}^{\prime}dx\mid$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}\mid
a_{11}u_{1h}\mid^{2}dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}\mid
a_{12}u_{2h}\mid^{2}dx+\frac{1}{2}\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}\mid
b_{1}u_{1h}^{\prime}\mid^{2}dx$ $\displaystyle\leq$
$\displaystyle\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}(\frac{\sigma}{4\mu^{2}})\mu^{2}\mid
u_{1h}\mid^{2}dx+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}(\frac{\sigma}{4\mu^{2}})\mu^{2}\mid
u_{2h}\mid^{2}dx+\frac{1}{2}\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}\mid
b_{1}u_{1h}^{\prime}\mid^{2}dx$ $\displaystyle=$
$\displaystyle\frac{\sigma}{4}(\|u_{1h}\|_{0}^{2}+\|u_{2h}\|_{0}^{2})+\frac{1}{2}\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(b_{1}u_{1h}^{\prime})^{2}dx$

and similarly we have

$\mid\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(a_{21}u_{1h}+a_{22}u_{2h})b_{2}u_{2h}^{\prime}dx\mid\leq\frac{\sigma}{4}(\|u_{1h}\|_{0}^{2}+\|u_{2h}\|_{0}^{2})+\frac{1}{2}\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(b_{2}u_{1h}^{\prime})^{2}dx.$

Combining the above two results we have the desired result. Hence $B_{h}$ is
coercive with respect to $|||.|||_{V_{h}}.$ ∎

Also $B_{h}$ is continuous in the discrete energy norm and $f_{h}$ is a
bounded linear functional on $V_{h}^{2}.$ By Lax-Milgram Theorem, we conclude
that the problem (2.4) has a unique solution.

######  Remark 2.3.

While deriving the corresponding difference scheme, we use the SDFEM with
lumping for the terms $(a_{11}u_{1}+a_{12}u_{2},v_{1})$ and
$(a_{21}u_{1}+a_{22}u_{2},v_{2}).$ That is $(a_{11}u_{1},v_{1})$ is replaced
by $\sum_{i=1}^{N-1}\bar{h}_{i}\widehat{a_{11,i}}u_{1,i}v_{1,i}$ where
$\widehat{a_{11,i}}=\frac{\bar{b}_{1}^{2}}{\beta_{1}^{2}}\|a_{11}\|_{\infty[x_{i},x_{i-1}]}.$

We choose $d=x_{N/2}$ and take
$f_{1}(d)=f_{1}(N/2)=\frac{f_{1}(\frac{N}{2}-1)+f_{1}(\frac{N}{2}+1)}{2},\,\quad$
$f_{2}(d)=f_{2}(N/2)=\frac{f_{2}(\frac{N}{2}-1)+f_{2}(\frac{N}{2}+1)}{2}.$
Then the corresponding difference scheme is

$\begin{split}L^{N}\bar{U}_{i}&:=\begin{cases}-\varepsilon[(\frac{U_{1,i+1}-U_{1,i}}{h_{i+1}}-\frac{U_{1,i}-U_{1,i-1}}{h_{i}})+(\frac{U_{2,i+1}-U_{2,i}}{h_{i+1}}-\frac{U_{2,i}-U_{2,i-1}}{h_{i}})]\\\
+\alpha_{1,i}(\frac{U_{1,i+1}-U_{1,i}}{h_{i+1}})+\alpha_{2,i}(\frac{U_{2,i+1}-U_{2,i}}{h_{i+1}})\\\
+\beta_{1,i}(\frac{U_{1,i}-U_{1,i-1}}{h_{i}})+\beta_{2,i}(\frac{U_{2,i}-U_{2,i-1}}{h_{i}})\\\
+\gamma_{1,i}U_{1,i}+\gamma_{2,i}U_{2,i}=f_{h}(\bar{\phi}_{i}),\end{cases}\\\
&U_{1,0}=U_{1,N}=U_{2,0}=U_{2,N}=0,\end{split}$ (2.5)

where $\bar{U}_{i}=(U_{1,i},U_{2,i}),\,$
$U_{1,i}=U_{1}(x_{i}),\,U_{2,i}=U_{2}(x_{i}),\,\bar{\phi}_{i}=(\phi_{i},\phi_{i}),\quad
i=1,2,...,N-1$ and

$\displaystyle\alpha_{1,i}$
$\displaystyle=h_{i+1}\int_{x_{i}}^{x_{i+1}}(b_{1}\phi_{i+1}^{\prime}\phi_{i}+\delta_{1,i+1}b_{1}^{2}\phi_{i+1}^{\prime}\phi_{i}^{\prime}+\delta_{1,i+1}b_{1}a_{11}\phi_{i+1}\phi_{i}^{\prime}+\delta_{2,i+1}b_{2}a_{21}\phi_{i+1}\phi_{i}^{\prime})dx$
$\displaystyle\beta_{1,i}$
$\displaystyle=-h_{i}\int_{x_{i-1}}^{x_{i}}(b_{1}\phi_{i-1}^{\prime}\phi_{i}+\delta_{1,i}b_{1}^{2}\phi_{i-1}^{\prime}\phi_{i}^{\prime}+\delta_{1,i}b_{1}a_{11}\phi_{i-1}\phi_{i}^{\prime}+\delta_{2,i}b_{2}a_{21}\phi_{i-1}\phi_{i}^{\prime})dx$
$\displaystyle\gamma_{1,i}$
$\displaystyle=\bar{h}_{i}(\widehat{a_{11}}+\widehat{a_{21}})(x_{i})+\int_{x_{i-1}}^{x_{i}}(\delta_{1,i}b_{1}a_{11}+\delta_{2,i}b_{2}a_{21})\phi_{i}^{\prime}dx+\int_{x_{i}}^{x_{i+1}}(\delta_{1,i+1}b_{1}a_{11}+\delta_{2,i+1}b_{2}a_{21})\phi_{i}^{\prime}dx$
$\displaystyle\alpha_{2,i}$
$\displaystyle=h_{i+1}\int_{x_{i}}^{x_{i+1}}(b_{2}\phi_{i+1}^{\prime}\phi_{i}+\delta_{2,i+1}b_{2}^{2}\phi_{i+1}^{\prime}\phi_{i}^{\prime}+\delta_{1,i+1}b_{1}a_{12}\phi_{i+1}\phi_{i}^{\prime}+\delta_{2,i+1}b_{2}a_{22}\phi_{i+1}\phi_{i}^{\prime})dx$
$\displaystyle\beta_{2,i}$
$\displaystyle=-h_{i}\int_{x_{i}}^{x_{i+1}}(b_{2}\phi_{i-1}^{\prime}\phi_{i}+\delta_{2,i}b_{2}^{2}\phi_{i-1}^{\prime}\phi_{i}^{\prime}+\delta_{1,i}b_{1}a_{12}\phi_{i-1}\phi_{i}^{\prime}+\delta_{2,i}b_{2}a_{22}\phi_{i-1}\phi_{i}^{\prime})dx$
$\displaystyle\gamma_{2,i}$
$\displaystyle=\bar{h}_{i}(\widehat{a_{12}}+\widehat{a_{22}})(x_{i})+\int_{x_{i-1}}^{x_{i}}(\delta_{1,i}b_{1}a_{12}+\delta_{2,i}b_{2}a_{22})\phi_{i}^{\prime}dx+\int_{x_{i}}^{x_{i+1}}(\delta_{1,i+1}b_{1}a_{12}+\delta_{2,i+1}b_{2}a_{22})\phi_{i}^{\prime}dx.$

######  Remark 2.4.

If the local mesh step is small enough, then it is possible to choose
$\delta_{k,i}=0,k=1,2.$ In other case, we shall choose $\delta_{k,i}$ from the
condition, $\alpha_{k,i}$ of the difference scheme (2.5) equal to zero. Thus
we have

$\displaystyle\delta_{1,i}=\begin{cases}0,\quad\quad\quad\quad\quad\quad\quad
h_{i}\leq\frac{2\varepsilon}{\parallel b_{1}\parallel_{\infty}},\\\
\frac{b_{1}h_{i}(2b_{2}^{2}+h_{i}b_{2}a_{22})-h_{i}^{2}b_{2}^{2}a_{21}}{(2b_{1}^{2}+h_{i}b_{1}a_{11})(2b_{2}^{2}+h_{i}b_{2}a_{22})-h_{i}^{2}b_{1}b_{2}a_{12}a_{21}},h_{i}>\frac{2\varepsilon}{\parallel
b_{1}\parallel_{\infty}}\end{cases}$

and also

$\displaystyle\delta_{2,i}=\begin{cases}0,\quad\quad\quad\quad\quad\quad\quad
h_{i}\leq\frac{2\varepsilon}{\parallel b_{2}\parallel_{\infty}},\\\
\frac{b_{2}h_{i}(2b_{1}^{2}+h_{i}b_{1}a_{11})-h_{i}^{2}b_{1}^{2}a_{12}}{(2b_{1}^{2}+h_{i}b_{1}a_{11})(2b_{2}^{2}+h_{i}b_{2}a_{22})-h_{i}^{2}b_{1}b_{2}a_{12}a_{21}},h_{i}>\frac{2\varepsilon}{\parallel
b_{2}\parallel_{\infty}}.\end{cases}$

We derive the following estimates of $\delta_{1,i}$ and $\delta_{2,i}$

$\displaystyle\delta_{k,i}\leq\begin{cases}CN^{-1}\quad\text{for}\quad
i=1,...,N/4\quad\text{and}\quad i=(N/2)+1,...,3N/4,\\\
0\quad\quad\quad\text{for}\quad i=(N/4)+1,...,N/2\quad\text{and}\quad
i=(3N/4)+1,...,N-1,\end{cases}$

where $k=1,2.$

The above system contains $N-1$ equations and has $2N-2$ unknowns. To solve
the system we split it into two algebraic systems as follows:
For $i=1,2,...,N-1$

$\displaystyle
P_{1}^{N}U_{1,i}^{*}:=\begin{cases}-\varepsilon(\frac{U_{1,i+1}^{*}-U_{1,i}^{*}}{h_{i+1}}-\frac{U_{1,i}^{*}-U_{1,i-1}^{*}}{h_{i}})+\alpha_{1,i}(\frac{U_{1,i+1}^{*}-U_{1,i}^{*}}{h_{i+1}})+\beta_{1,i}(\frac{U_{1,i}^{*}-U_{1,i-1}^{*}}{h_{i}})\\\
+\gamma_{1,i}U_{1,i}^{*}=\int_{x_{i-1}}^{x_{i+1}}f_{1}\phi_{i}+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}f_{1}b_{1}\phi_{i}^{\prime},\quad
U_{1,0}^{*}=U_{1,N}^{*}=0,\end{cases}$ (2.6) $\displaystyle
P_{2}^{N}U_{2,i}^{*}:=\begin{cases}-\varepsilon(\frac{U_{2,i+1}^{*}-U_{2,i}^{*}}{h_{i+1}}-\frac{U_{2,i}^{*}-U_{2,i-1}^{*}}{h_{i}})+\alpha_{2,i}(\frac{U_{2,i+1}^{*}-U_{2,i}^{*}}{h_{i+1}})+\beta_{2,i}(\frac{U_{2,i}^{*}-U_{2,i-1}^{*}}{h_{i}})\\\
+\gamma_{2,i}U_{2,i}^{*}=\int_{x_{i-1}}^{x_{i+1}}f_{2}\phi_{i}+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}f_{2}b_{2}\phi_{i}^{\prime},\quad
U_{2,0}^{*}=U_{2,N}^{*}=0.\end{cases}$ (2.7)

The above system (2.6) corresponds to the differential equation

$P_{1}^{*}u_{1}^{*}:=-\varepsilon
u_{1}^{*^{\prime\prime}}+b_{1}(x)u_{1}^{*^{\prime}}+(a_{11}(x)+a_{21}(x))u_{1}^{*}=f_{1}(x),\quad
x\in(\Omega^{-}\cup\Omega^{+}),$

subject to boundary conditions $u_{1}^{*}(0)=u_{1}^{*}(1)=0.$ This boundary
value problem has a unique solution [6]. Using the inverse monotone property
of the matrix, one can establish the numerical stability of the system (2.6).
Similarly we can deal with second equation (2.7). If $U_{1,i}^{*}$ and
$U_{2,i}^{*}$ are solutions of (2.6) and (2.7) respectively then
$(U_{1,i}^{*},U_{2,i}^{*})$ is a solution of (2.5). By uniqueness, this is the
only possible solution. Therefore, it is enough to solve (2.6) and (2.7).

##  3\. Error analysis

The convergence analysis of the numerical scheme starts at the triangle
inequality

$|||\bar{u}-\bar{u}_{h}|||_{V_{h}}\leq|||\bar{u}-\bar{u}^{I}|||_{V_{h}}+|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}},$
(3.1)

where $\bar{u}^{I}$ denotes the piecewise linear interpolant to $\bar{u}$ on
$\Omega$.
Now we estimate the second term of equation (3.1).

######  Lemma 3.1.

The following estimate holds true

$|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}}\leq C\|\bar{u}^{I}-\bar{u}\|_{0}.$

###### Proof.

Because of the Galerkin orthogonality relation between $\bar{u}$ and
$\bar{u}_{h}$, we have

$B_{h}(\bar{u}_{h}-\bar{u},\bar{u}^{I}-\bar{u}_{h})=0.$

Then from the coercive property (2.2) of $B_{h}(.,.),$ we have

$\displaystyle|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}}^{2}$ $\displaystyle\leq$
$\displaystyle 2B_{h}(\bar{u}^{I}-\bar{u}_{h},\bar{u}^{I}-\bar{u}_{h})$
$\displaystyle=$ $\displaystyle
2B_{h}(\bar{u}^{I}-\bar{u},\bar{u}^{I}-\bar{u}_{h})$ $\displaystyle=$
$\displaystyle
2[(b_{1}(u_{1}^{I}-u_{1})^{\prime},u_{1}^{I}-u_{1h})+(a_{11}(u_{1}^{I}-u_{1}),u_{1}^{I}-u_{1h})+(a_{12}(u_{2}^{I}-u_{2}),u_{1}^{I}-u_{1h})$
$\displaystyle+(b_{2}(u_{2}^{I}-u_{2})^{\prime},u_{2}^{I}-u_{2h})+(a_{21}(u_{1}^{I}-u_{1}),u_{2}^{I}-u_{2h})+(a_{22}(u_{2}^{I}-u_{2}),u_{2}^{I}-u_{2h})$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}(-\varepsilon(u_{1}^{I}-u_{1})^{\prime\prime}+b_{1}(u_{1}^{I}-u_{1})^{\prime}+a_{11}(u_{1}^{I}-u_{1})$
$\displaystyle+a_{12}(x)(u_{2}^{I}-u_{2}))b_{1}(u_{1}^{I}-u_{1h})^{\prime}dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}(-\varepsilon(u_{2}^{I}-u_{2})^{\prime\prime}+b_{2}(u_{2}^{I}-u_{2})^{\prime}+a_{21}(u_{1}^{I}-u_{1})$
$\displaystyle+a_{22}(u_{2}^{I}-u_{2}))b_{2}(u_{2}^{I}-u_{2h})^{\prime}dx].$

That is,

$\displaystyle|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}}^{2}$ $\displaystyle\leq$
$\displaystyle
C\int_{0}^{1}(u_{1}^{I}-u_{1})[(u_{1}^{I}-u_{1h})+(u_{2}^{I}-u_{2h})]+C\int_{0}^{1}(u_{2}^{I}-u_{2})[(u_{1}^{I}-u_{1h})+(u_{2}^{I}-u_{2h})]$
$\displaystyle\leq$ $\displaystyle
C\int_{0}^{1}[(u_{1}^{I}-u_{1})+(u_{2}^{I}-u_{2})][(u_{1}^{I}-u_{1h})+(u_{2}^{I}-u_{2h})].$

Therefore we have

$\displaystyle|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}}^{2}$ $\displaystyle\leq$
$\displaystyle
C\|\bar{u}^{I}-\bar{u}\|_{0}\quad\|\bar{u}^{I}-\bar{u}_{h}\|_{0}$
$\displaystyle|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}}^{2}$ $\displaystyle\leq$
$\displaystyle
C\|\bar{u}^{I}-\bar{u}\|_{0}\quad|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}}$
$\displaystyle|||\bar{u}^{I}-\bar{u}_{h}|||_{V_{h}}$ $\displaystyle\leq$
$\displaystyle C\|\bar{u}^{I}-\bar{u}\|_{0}.$

∎

###  3.1. Error analysis on Shishkin and Bakhvalov-Shishkin meshes

For the discretization described above we shall use a mesh of the general type
introduced in [9], but here adapted for the layers at $x=d.$ Let $N>4$ be a
positive even integer and

$\sigma_{1}=\min\\{\frac{d}{2},{\frac{\varepsilon}{\beta}}\tau_{0}\ln
N\\},\quad\sigma_{2}=\min\\{\frac{1-d}{2},{\frac{\varepsilon}{\beta}}\tau_{0}\ln
N\\},\quad\quad\tau_{0}\geq 2.$

Our mesh will be equidistant on $\bar{\Omega}_{S}$, where

$\Omega_{S}=(0,d-\sigma_{1})\cup(d,1-\sigma_{2})$

and graded on $\bar{\Omega}_{0}$ where

$\Omega_{0}=(d-\sigma_{1},d)\cup(1-\sigma_{2},1).$

First we shall assume
$\sigma_{1}=\sigma_{2}=\frac{\tau_{0}\varepsilon}{\beta}\ln N$ as otherwise
$N^{-1}$ is exponentially small compared to $\varepsilon.$ We choose the
transition points to be

$x_{N/4}=d-\sigma_{1},\quad x_{N/2}=d,\quad x_{3N/4}=1-\sigma_{2}.$

Because of the specific layers, here we have to use two mesh generating
functions $\varphi_{1}$ and $\varphi_{2}$ which are both continuous and
piecewise continuously differentiable and monotonically decreasing functions
and

$\displaystyle\varphi_{1}(1/4)=\ln N,$ $\displaystyle\varphi_{1}(1/2)=0$
$\displaystyle\varphi_{2}(3/4)=\ln N,$ $\displaystyle\varphi_{2}(1)=0.$

The mesh points are

$x_{i}=\begin{cases}\frac{4i}{N}(d-\sigma_{1}),\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
i=0,...,N/4\\\
d-\frac{\tau_{0}}{\beta}\varepsilon\varphi_{1}(t_{i}),\quad\quad\quad\quad\quad\quad\quad\quad\quad
i=N/4+1,...,N/2\\\ d+\frac{4}{N}(1-d-\sigma_{2})(i-N/2),\quad\quad\quad
i=N/2+1,...,3N/4\\\
1-\frac{\tau_{0}}{\beta}\varepsilon\varphi_{2}(t_{i}),\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
i=3N/4+1,...,N,\end{cases}$

where $t_{i}=i/N$. We define new functions $\psi_{1}$ and $\psi_{2}$ by

$\varphi_{i}=-\ln\psi_{i},\quad i=1,2.$

There are several mesh-characterizing functions $\psi$ in the literature, but
we shall use only those which correspond to Shishkin mesh and Bakhvalov-
Shishkin mesh with the following properties

$\displaystyle\max|\psi^{\prime}|$ $\displaystyle=$ $\displaystyle C\ln
N\quad\text{for Shishkin meshes}$ $\displaystyle\max|\psi^{\prime}|$
$\displaystyle=$ $\displaystyle C\quad\text{for Bakhvalov-Shishkin meshes}$

$\bullet$ Shishkin mesh

$\psi_{1}(t)=e^{-2(1-2t)lnN},\quad\,\ \psi_{2}(t)=e^{-4(1-t)lnN},$

$\bullet$ Bakhvalov-Shishkin mesh

$\psi_{1}(t)=1-2(1-N^{-1})(1-2t),\quad\psi_{2}(t)=1-4(1-N^{-1})(1-t).$

The set of interior mesh points is denoted by
$\Omega_{\varepsilon}^{N}=\bar{\Omega}_{\varepsilon}^{N}\setminus\\{x_{N/2}\\}$.
Also, for the both meshes, on the coarse part $\Omega_{S}$ we have

$h_{i}\leq CN^{-1}.$

It is well known that on the layer part of the Shishkin mesh [6]

$h_{i}\leq C\varepsilon N^{-1}\ln N$

and of the Bakhvalov-Shishkin mesh we have

$h_{i}\leq\begin{cases}\frac{\tau_{0}}{\beta}\varepsilon
N^{-1}\max{\mid\psi_{1}^{\prime}\mid}\exp{(\frac{\beta}{\tau_{0}\varepsilon}(d-x_{i-1}))},\quad
i=N/4+1,...,N/2,\\\ \frac{\tau_{0}}{\beta}\varepsilon
N^{-1}\max{\mid\psi_{2}^{\prime}\mid}\exp{(\frac{\beta}{\tau_{0}\varepsilon}(1-x_{i-1}))},\quad
i=3N/4+1,...,N\end{cases}$

and

$\frac{h_{i}}{\varepsilon}\leq CN^{-1}\max{\mid\varphi^{\prime}\mid}\leq C.$

##  4\. Interpolation Error

Initially we consider the interpolation error in the maximum norm. Let $f\in
C^{2}[x_{i-1},x_{i}]$ be arbitrary and $f^{I}$ a piecewise linear interpolant
to $f$ on $\Omega$. Then from the classical theory, we have

$|(f^{I}-f)(x)|\leq
2\int_{x_{i-1}}^{x_{i}}|f^{\prime\prime}(t)|(t-x_{i-1})dt.$

Now we compute the interpolation error for the first component $u_{1}$.

######  Lemma 4.1.

For the Shishkin mesh we have

$|u_{i}(x)-u_{i}^{I}(x)|\leq\begin{cases}CN^{-2}\ln^{2}N,x\in\Omega_{0}\\\
CN^{-2},x\in\Omega_{S}\end{cases}$

and for the Bakhavalov-Shishkin mesh it holds

$|u_{i}(x)-u_{i}^{I}(x)|\leq CN^{-2},x\in\Omega^{-}\cup\Omega^{+},\quad
i=1,2.$

###### Proof.

We now give a proof for the case $i=1$ for the Shishkin mesh. To prove the
estimates we use the decomposition of solution as smooth and layer components
and triangle inequality

$|(u_{1}-u_{1}^{I})(x)|\leq|(v_{1}-v_{1}^{I})(x)|+|(w_{11}-w_{11}^{I})(x)|+|(w_{21}-w_{21}^{I})(x)|.$
(4.1)

On Shishkin meshes, let $x\in[x_{i-1},x_{i}]\subset\Omega^{-}\cap\Omega_{S}.$
Then the first term of (4.1) will be

$\displaystyle|(v_{1}-v_{1}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
2\int_{x_{i-1}}^{x_{i}}|v_{1}^{\prime\prime}(t)|(t-x_{i-1})dt$
$\displaystyle\leq$ $\displaystyle 2C\int_{x_{i-1}}^{x_{i}}(t-x_{i-1})dt$
$\displaystyle\leq$ $\displaystyle 2C\frac{h_{i}^{2}}{2}$
$\displaystyle|(v_{1}-v_{1}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
CN^{-2}.$

Again the second term of (4.1) will be

$\displaystyle|(w_{11}-w_{11}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
2\|w_{11}(x)\|_{L_{\infty}[x_{i-1},x_{i}]}$ $\displaystyle\leq$ $\displaystyle
C\max_{i}e^{\frac{-\beta(1-x_{i})}{\varepsilon}}$
$\displaystyle|(w_{11}-w_{11}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
CN^{-\tau_{0}}.$

To compute the last term of (4.1), we have

$\displaystyle|(w_{21}-w_{21}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
2\|w_{21}(x)\|_{L_{\infty}[x_{i-1},x_{i}]}$ $\displaystyle\leq$ $\displaystyle
C\varepsilon\max_{i}e^{\frac{-\beta(d-x_{i})}{\varepsilon}}$
$\displaystyle\leq$ $\displaystyle
CN^{-1}\max_{i}e^{\frac{-\beta(d-x_{i})}{\varepsilon}}$
$\displaystyle|(w_{21}-w_{21}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
CN^{-1-\tau_{0}}.$

Now let $x\in[x_{i-1},x_{i}]\subset\Omega^{-}\cap\Omega_{0}$ we have

$\displaystyle|(v_{1}-v_{1}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
2\int_{x_{i-1}}^{x_{i}}|v_{1}^{\prime\prime}(t)|(t-x_{i-1})dt$
$\displaystyle\leq$ $\displaystyle 2C\int_{x_{i-1}}^{x_{i}}(t-x_{i-1})dt$
$\displaystyle\leq$ $\displaystyle C\frac{h_{i}^{2}}{2}$ $\displaystyle\leq$
$\displaystyle C(\varepsilon N^{-1}\ln N)^{2}$

and also the second term on $\Omega_{0}$ will be

$\displaystyle|(w_{11}-w_{11}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
2\|w_{11}(x)\|_{L_{\infty}[x_{i-1},x_{i}]}$ $\displaystyle\leq$ $\displaystyle
2\max_{i}e^{\frac{-\beta(1-x_{i})}{\varepsilon}}$
$\displaystyle|(w_{11}-w_{11}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
CN^{-\tau_{0}}.$

The last term on $\Omega_{0}$ will be

$\displaystyle|(w_{21}-w_{21}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
2\|w_{21}(x)\|_{L_{\infty}[x_{i-1},x_{i}]}$ $\displaystyle\leq$ $\displaystyle
C\varepsilon\max_{i}e^{\frac{-\beta(d-x_{i})}{\varepsilon}}$
$\displaystyle\leq$ $\displaystyle
CN^{-1}\max_{i}e^{\frac{-\beta(d-x_{i})}{\varepsilon}}.$
$\displaystyle|(w_{21}-w_{21}^{I})(x)|$ $\displaystyle\leq$ $\displaystyle
CN^{-1-\tau_{0}}.$

Similarly we will also obtain the same estimate on $x\in\Omega^{+}.$ From
equation (4.1), hence the result.
On Bakhavalov-Shishkin mesh, we follow the above similar procedure to obtain
the result. ∎

Now we consider the interpolation error of $\bar{u}$ in $L_{2}$-norm

$\|\bar{u}-\bar{u}^{I}\|_{0}=[(\int_{0}^{1}|u_{1}-u_{1}^{I}|^{2}dx)+(\int_{0}^{1}|u_{2}-u_{2}^{I}|^{2}dx)]^{1/2}.$
(4.2)

######  Lemma 4.2.

For Shishkin mesh, the interpolation error of $\bar{u}$ in $L_{2}$-norm is

$\|\bar{u}-\bar{u}^{I}\|_{0}\leq CN^{-5/2}\ln^{5/2}N.$

###### Proof.

Consider the first component of equation (4.2)

$\displaystyle\int_{0}^{1}|u_{1}-u_{1}^{I}|^{2}dx$ $\displaystyle\leq$
$\displaystyle\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}|u_{1}-u_{1}^{I}|^{2}dx$
$\displaystyle\leq$ $\displaystyle\sum_{i=1}^{N}(CN^{-2}\ln^{2}N)^{2}h_{i}$
$\displaystyle\leq$ $\displaystyle C_{1}(CN^{-2}\ln^{2}N)^{2}(C\varepsilon\ln
N)$ $\displaystyle\|u_{1}-u_{1}^{I}\|_{0}$ $\displaystyle\leq$ $\displaystyle
CN^{-5/2}\ln^{5/2}N.$

Similarly one can easily prove

$\|u_{2}-u_{2}^{I}\|_{0}\leq CN^{-5/2}\ln^{5/2}N.$

From (4.2), we have an estimate of $\bar{u}-\bar{u}^{I}$ in $L_{2}-$ norm

$\|\bar{u}-\bar{u}^{I}\|_{0}\leq CN^{-5/2}\ln^{5/2}N.$

∎

######  Lemma 4.3.

Let $\bar{u}$ and $\bar{u}^{I}$ be solution of (1.1-1.3) and linear
interpolant of $\bar{u}$ respectively. Then we have

$|||\bar{u}-\bar{u}^{I}|||_{V_{h}}\leq CN^{-1}\ln^{3/2}N,\text{for Shishkin
meshes}$

###### Proof.

Since

$\int_{0}^{1}((u_{1}-u_{1}^{I})^{\prime}(x))^{2}dx=-\int_{0}^{1}(u_{1}-u_{1}^{I})(x)u_{1}^{\prime\prime}(x)dx$

therefore, by Lemma 4.1 we conclude that

$\displaystyle\int_{0}^{1}|(u_{1}-u_{1}^{I})^{\prime}(x)|^{2}dx$
$\displaystyle\leq$ $\displaystyle
C\max_{x_{i}\in\Omega_{\varepsilon}^{N}}|(u_{1}-u_{1}^{I})(x_{i})|\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}u_{1}^{{}^{\prime\prime}}(x)dx$
$\displaystyle\leq$ $\displaystyle
CN^{-2}\ln^{2}N\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}(v_{1}^{{}^{\prime\prime}}(x)+w_{11}^{{}^{\prime\prime}}(x)+w_{21}^{{}^{\prime\prime}}(x))dx.$

then for the regular part of the solution we have

$|\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}v_{1}^{{}^{\prime\prime}}(x)|\leq
C(\varepsilon\ln N+1)$

and for the singular part

$\displaystyle|\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}w_{11}^{{}^{\prime\prime}}(x)|$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-2}\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}e^{-\frac{\beta(1-x)}{\varepsilon}}dx$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-1}[\sum_{i=\frac{N}{4}+1}^{\frac{N}{2}}[e^{-\frac{\beta(1-x)}{\varepsilon}}]_{x_{i-1}}^{x_{i}}+\sum_{i=\frac{3N}{4}+1}^{N}[e^{-\frac{\beta(1-x)}{\varepsilon}}]_{x_{i-1}}^{x_{i}}]+C\varepsilon^{-1}N^{1-\tau_{0}}$
$\displaystyle\leq$ $\displaystyle C\varepsilon^{-1}\ln N.$

and

$\displaystyle|\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}w_{21}^{{}^{\prime\prime}}(x)|$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-1}[\sum_{i=1}^{\frac{N}{2}}\int_{x_{i-1}}^{x_{i}}e^{-\frac{\beta(d-x)}{\varepsilon}}dx+\sum_{i=\frac{N}{2}+1}^{N}\int_{x_{i-1}}^{x_{i}}e^{-\frac{\beta(1-x)}{\varepsilon}}dx]$
$\displaystyle\leq$ $\displaystyle CN^{1-\tau_{0}}$

Using the assumption $\tau_{0}\geq 2$ and the above estimates we have

$\int_{0}^{1}|(u_{1}-u_{1}^{I})^{\prime}(x)|^{2}dx\leq
CN^{-2}\ln^{2}N(\varepsilon\ln N+1+\varepsilon^{-1}\ln N+N^{-1})$

We also have similar result for $u_{2}$

$\int_{0}^{1}|(u_{2}-u_{2}^{I})^{\prime}(x)|^{2}dx\leq
C\varepsilon^{-1}N^{-2}\ln^{3}N.$

Now we combine the above results together

$|u_{1}-u_{1}^{I}|_{1}^{2}+|u_{2}-u_{2}^{I}|_{1}^{2}\leq
C\varepsilon^{-1}N^{-2}\ln^{3}N.$

Here we have to compute the interpolation error of $\bar{u}$ in energy norm,
that is, $|||\bar{u}-\bar{u}^{I}|||_{V_{h}}.$
We have

$\displaystyle|||\bar{u}-\bar{u}^{I}|||_{V_{h}}$
$\displaystyle=[\varepsilon(|u_{1}-u_{1}^{I}|_{1}^{2}+|u_{2}-u_{2}^{I}|_{1}^{2})+\sigma(\|u_{1}-u_{1}^{I}\|_{0}^{2}+\|u_{2}-u_{2}^{I}\|_{0}^{2})$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}b_{1}^{2}(x_{i})((u_{1}-u_{1}^{I})^{\prime}(x))^{2}dx$
$\displaystyle+\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}b_{2}^{2}(x_{i})((u_{2}-u_{2}^{I})^{\prime}(x))^{2}dx]^{1/2}.$

Now we have to estimate the following terms

$\displaystyle|\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}b_{1}^{2}(x_{i})((u_{1}-u_{1}^{I})^{\prime}(x))^{2}dx|$
$\displaystyle\leq
C\mid\delta_{1,i}\mid\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\mid(u_{1}-u_{1}^{I})^{\prime}(x)^{2}\mid
dx$ $\displaystyle\leq CN^{-1}(\varepsilon^{-1}N^{-2}\ln^{3}N)$
$\displaystyle|\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{1,i}b_{1}^{2}(x_{i})((u_{1}-u_{1}^{I})^{\prime}(x))^{2}dx|$
$\displaystyle\leq CN^{-2}\ln^{3}N$

and also

$\displaystyle|\sum_{i=1}^{N}\int_{x_{i-1}}^{x_{i}}\delta_{2,i}b_{2}^{2}(x_{i})((u_{2}-u_{2}^{I})^{\prime}(x))^{2}dx|\leq
CN^{-2}\ln^{3}N.$

Substituting these estimates, we have

$\displaystyle|||\bar{u}-\bar{u}^{I}|||_{V_{h}}$
$\displaystyle\leq[\varepsilon(C\varepsilon^{-1}N^{-2}\ln^{3}N)+\sigma(CN^{-5/2}\ln^{5/2}N)^{2}+CN^{-3}\ln^{4}N+CN^{-2}\ln^{3}N]^{1/2}$
$\displaystyle\leq[\varepsilon(C\varepsilon^{-1}N^{-2}\ln^{3}N)+\sigma(CN^{-5/2}\ln^{5/2}N)^{2}+CN^{-2}\ln^{3}N]^{1/2}$
$\displaystyle\leq CN^{-1}\ln^{3/2}N[1+N^{-3}\ln^{2}N+1]^{1/2}$
$\displaystyle|||\bar{u}-\bar{u}^{I}|||_{V_{h}}$ $\displaystyle\leq
CN^{-1}\ln^{3/2}N.$

∎

##  5\. Error Estimate

Now we state the main theorem of this paper.

######  Theorem 5.1.

Let $\bar{u}$ and $\bar{u}_{h}$ be solution of (1.1-1.3) and (2.4)
respectively. Then we have

$|||\bar{u}-\bar{u}_{h}|||_{V_{h}}\leq\begin{cases}CN^{-1}\ln^{3/2}N,\quad\text{for
Shishkin mesh,}\\\ CN^{-1},\quad\text{for Bakhvalov-Shishkin
mesh}.\end{cases}$

###### Proof.

From the inequality (3.1), Lemmas (3.1), (4.2) and (4.3), for Shishkin meshes
we have

$\displaystyle|||\bar{u}-\bar{u}_{h}|||_{V_{h}}$ $\displaystyle\leq$
$\displaystyle CN^{-1}\ln^{3/2}N+CN^{-\frac{5}{2}}\ln^{\frac{5}{2}}N$
$\displaystyle\leq$ $\displaystyle CN^{-1}\ln^{3/2}N.$

Similarly we prove the error estimates for Bakhvalov-Shishkin meshes. ∎

##  6\. Numerical Experiments

In this section we experimentally verify our theoretical results proved in the
previous section.

######  Example 6.1.

Consider the BVP

$\displaystyle-\varepsilon
u_{1}^{{}^{\prime\prime}}(x)+u_{1}^{\prime}(x)+2u_{1}(x)-u_{2}(x)=f_{1}(x),\quad
x\in\Omega^{-}\cup\Omega^{+},$ (6.1) $\displaystyle-\varepsilon
u_{2}^{{}^{\prime\prime}}(x)+u_{2}^{\prime}(x)-u_{1}(x)+2u_{2}(x)=f_{2}(x),\quad
x\in\Omega^{-}\cup\Omega^{+},$ (6.2) $u_{1}(0)=0,\quad u_{1}(1)=0,\quad
u_{2}(0)=0,\quad u_{2}(1)=0,$ (6.3)

where

$f_{1}(x)=\begin{cases}1,\quad 0\leq x\leq 0.5,\\\ -0.8,\quad 0.5\leq x\leq
1\end{cases}$

and

$f_{2}(x)=\begin{cases}-2.0,\quad 0\leq x\leq 0.5,\\\ 1.8,\quad 0.5\leq x\leq
1\end{cases}$

For our tests, we take $\varepsilon=2^{-18}$, which is sufficiently small to
bring out the singularly perturbed nature of the problem. Now we define a
maximum norm of $\bar{u}_{h}$ as

$\parallel\bar{u}_{h}\parallel_{\infty}=\max\\{\max_{1\leq i\leq N-1}\\{\mid
u_{1h}(x_{i})\mid\\},\max_{1\leq i\leq N-1}\\{\mid u_{2h}(x_{i})\mid\\}\\}$

We measure the accuracy in various norms and the rates of convergence $r^{N}$
are computed using the following formula:

$r^{N}=log_{2}(\dfrac{E^{N}}{E^{2N}}),$

where

$E^{N}=\begin{cases}\parallel\bar{u}_{h}-\bar{u}_{2h}^{I}\parallel_{\infty},\quad\text{for
maximum norm,}\\\
\parallel\bar{u}_{h}-\bar{u}_{2h}^{I}\parallel_{0},\quad\text{for}(L_{2}(\Omega))^{2}-\text{norm,}\\\
|||\bar{u}_{h}-\bar{u}_{2h}^{I}|||_{V_{h}},\quad\text{ for discrete energy
norm},\end{cases}$

and $\bar{u}_{h}^{I}$ denotes the piecewise linear interpolant of $\bar{U}.$
In Tables 1 and 2, we present values of $E^{N},r^{N}$ for the solution of the
BVP (6.1)-(6.3) for Shishkin and Bakhvalov-Shishkin meshes respectively. The
Figures 1 and 2 depict the numerical solution of the BVP (6.1)-(6.3) for
Shishkin mesh. We compare the values of $E^{N},r^{N}$ for the solution of the
same BVP (6.1)-(6.3) for Shishkin mesh using the standard upwind scheme
adopted [13]. From the tables, we infer that the order of convergence is
higher in the cases of maximum norm and $L_{2}-$ norm when compared with
discrete energy norm as defined earlier. Therefore the present method may
yield better results.
The numerical results are clear illustrations of the convergence estimates
derived in the present paper for both the type of meshes.

######  Remark 6.2.

It may be observed that the value of $\tau_{0}$ is taken as $\tau_{0}\geq 2.$
From the above experimental results this condition seems to be essential.
Infact, it is found that if one takes the value $\tau_{0}<2$ the order of
convergence may not be $2.$

Table 1. Values of $E^{N}$ and $r^{N}$ for the solution of the BVP (6.1) - (6.3) in different norms for Shishkin mesh. N | $\parallel\bar{u}_{h}-\bar{u}_{h}^{I}\parallel_{\infty}$ | $\parallel\bar{u}_{h}-\bar{u}_{h}^{I}\parallel_{0}$ | $|||\bar{u}_{h}-\bar{u}_{h}^{I}|||_{V_{h}}$
---|---|---|---
| $E^{N}$ | $r^{N}$ | $E^{N}$ | $r^{N}$ | $E^{N}$ | $r^{N}$
$32$ | 2.3693e-01 | 1.4253 | 1.0785e-02 | 1.1113 | 2.7108e-01 | 0.8742
$64$ | 8.8222e-02 | 1.0592 | 4.9921e-03 | 1.0447 | 1.4788e-01 | 0.6716
$128$ | 4.2337e-02 | 0.9939 | 2.4199e-03 | 1.0176 | 9.2838e-02 | 0.5973
$256$ | 2.1258e-02 | 0.9986 | 1.1953e-03 | 1.0061 | 6.1517e-02 | 0.5625
$512$ | 1.0639e-02 | 1.0030 | 5.9513e-04 | 1.0011 | 4.1654e-02 | 0.5620
$1024$ | 5.3085e-03 | 1.0085 | 2.9734e-04 | 0.9990 | 2.8213e-02 | 0.5522
$2048$ | 2.6387e-03 | - | 1.4877e-04 | - | 1.9240e-02 | -
Table 2. Values of $E^{N}$ and $r^{N}$ for the solution of the BVP (6.1) - (6.3) in different norms for Bakhvalov-Shishkin mesh. N | $\parallel\bar{u}_{h}-\bar{u}_{h}^{I}\parallel_{\infty}$ | $\parallel\bar{u}_{h}-\bar{u}_{h}^{I}\parallel_{0}$ | $|||\bar{u}_{h}-\bar{u}_{h}^{I}|||_{V_{h}}$
---|---|---|---
| $E^{N}$ | $r^{N}$ | $E^{N}$ | $r^{N}$ | $E^{N}$ | $r^{N}$
$32$ | 1.6550-01 | 0.9811 | 1.1047e-02 | 0.8554 | 2.7386e-01 | 0.5465
$64$ | 8.3838e-02 | 0.9865 | 4.9717e-03 | 0.9120 | 1.8750e-01 | 0.5304
$128$ | 4.2313e-02 | 0.9945 | 2.3671e-03 | 0.9535 | 1.2981e-01 | 0.5194
$256$ | 2.1236e-02 | 1.0001 | 1.1551e-03 | 0.9769 | 9.0558e-02 | 0.5157
$512$ | 1.0617e-02 | 1.0059 | 5.7064e-04 | 0.9874 | 6.3341e-02 | 0.5192
$1024$ | 5.2870e-03 | 1.0141 | 2.8361e-04 | 0.9940 | 4.4194e-02 | 0.5322
$2048$ | 2.6177e-03 | - | 1.4139e-04 | - | 3.0560e-02 | -

Figure 1. Graphs of the numerical solution of the first component $u_{1h}$ of
the BVP (6.1)-(6.3) for various values of $\varepsilon$ with $N=512.$

Figure 2. Graphs of the numerical solution of the second component $u_{2h}$ of
the BVP (6.1)-(6.3) for various values of $\varepsilon$ with $N=512.$

Figure 3. Plots of order of convergence for Example 6.1 and
$\varepsilon=2^{-18}$ in various norms.

## References

  * [1] P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O’Riordan, G. I. Shishkin, _Robust computational techniques for boundary layers_ , Chapman Hall/ CRC, Boca Raton, 2000.
  * [2] H-G. Roos, M. Stynes, L. Tobiska, _Numerical methods for singularly perturbed differential equations_ , Volume 24 of Springer series in Computational Mathematics, Springer-Verlag, Berlin, 1996.
  * [3] E.P. Doolan, J.J.H.Miller, W. H. A. Schilders, _Uniform numerical methods for problems with initial and boundary layers_ , Boole, Dublin, 1980.
  * [4] A. H. Nayfeh, _Introduction to Perturbation Methods_ , Wiley, New York, 1981.
  * [5] R. E. O’Malley, _Singular perturbation methods for ordinary differential equations_ , Springer, New York, 1990.
  * [6] H-G. Roos, Helena Zarin, _The streamline-diffusion method for a convection-diffusion problem with a point source,_ J. Comp. Appl. Math.,Vol.10, No.4(2002) 275-289.
  * [7] T. Lin$\ss{}$, N. Madden, _A finite element analysis of a coupled system of singularly perturbed reaction-diffusion equations_ , Appl. Math. Comput. 148(3)(2004)869-880.
  * [8] N. Madden, M. Stynes, _A uniformly convergent numerical method for a coupled system of two singularly perturbed linear reaction-diffusion problems_ , IMA. J. Numer. Anal.23(4)(2003)627-644.
  * [9] H-G. Roos, T. Lin$\ss{}$, _Sufficient conditions for uniform convergence on layer-adapted grids_ , Computing 63 (1999) 27-45.
  * [10] T. Lin$\ss{}$, N. Madden, _An improved error estimate for a system of coupled singularly perturbed reaction-diffusion equations_ , Comput. Methods Appl. Math. 3(2003)417-423.
  * [11] Zhongdi Cen, _Parameter-uniform finite difference scheme for a system of coupled singularly perturbed convection-diffusion equations_ , International Journal of Computer Mathematics, 82 (2005) 177-192.
  * [12] S. Bellew, E. O’Riordan, _A Parameter robust numerical method for a system of two singularly perturbed convection-diffusion equations_ , Journal of Applied Numerical Mathematrics, 51(2004) 171-186.
  * [13] A. Tamilselvan, N. Ramanujam, _A numerical method for singularly perturbed system of second order ordinary differential equations of convection-diffusion type with a discontinuous source term_ , Submitted for publication.