Further author information: (Send correspondence to Benjamin L. Gerard)
Benjamin L. Gerard: E-mail<EMAIL_ADDRESS>
# Various Wavefront Sensing and Control Developments on the Santa Cruz Extreme
AO Laboratory (SEAL) Testbed

Benjamin L. Gerard  UC Observatories  University of California Santa Cruz
Javier Perez-Soto  University of California Santa Cruz  Vincent Chambouleyron
Marseille Astrophysics Laboratory  Maaike A.M. van Kooten  UC Observatories
University of California Santa Cruz  Daren Dillon  UC Observatories  Sylvain
Cetre  Durahm University  Wakea Consulting  Rebecca Jensen-Clem  UC
Observatories  University of California Santa Cruz  Qiang Fu  King Abdullah
University of Science and Technology  Hadi Amata  King Abdullah University of
Science and Technology  Wolfgang Heidrich  King Abdullah University of Science
and Technology

###### Abstract

Ground-based high contrast imaging (HCI) and extreme adaptive optics (AO)
technologies have advanced to the point of enabling direct detections of gas-
giant exoplanets orbiting beyond the snow lines around nearby young star
systems. However, leftover wavefront errors using current HCI and AO
technologies, realized as “speckles” in the coronagraphic science image, still
limit HCI instrument sensitivities to detecting and characterizing lower-mass,
closer-in, and/or older/colder exoplanetary systems. Improving the performance
of AO wavefront sensors (WFSs) and control techniques is critical to improving
such HCI instrument sensitivity. Here we present three different ongoing
wavefront sensing and control project developments on the Santa cruz Extreme
AO Laboratory (SEAL) testbed: (1) “multi-WFS single congugate AO (SCAO)” using
the Fast Atmospheric Self-coherent camera (SCC) Technique (FAST) and a Shack
Hartmann WFS, (2) pupil chopping for focal plane wavefront sensing, first with
an external amplitude modulator and then with the DM as a phase-only
modulator, and (3) a laboratory demonstration of enhanced linearity with the
non-modulated bright Pyramid WFS (PWFS) compared to the regular PWFS. All
three topics share a common theme of multi-WFS SCAO and/or second stage AO,
presenting opportunities and applications to further investigate these
techniques in the future.

###### keywords:

Adaptive Optics, Wavefront Sensing, Focal Plane Wavefront Sensing, Pupil Plane
Wavefront Sensing

##  1 INTRODUCTION

Adaptive Optics (AO) technology development is entering a new era,
particularly in application to ground-based exoplanet and high contrast
imaging (HCI) science. Current state-of-the-art extreme AO (ExAO)-fed HCI
instruments enable near-infrared imaging of self-luminous exoplanets
$\gtrsim$a few Jupiter masses that are on $\gtrsim$Saturn-scale orbital
separations, and around nearby ($\lesssim$100 pc) young ($\lesssim$100 Myr)
host stars[1]. Future ExAO and HCI technology advancements are needed to
enable self-luminous gas giant exoplanet imaging down to Jupiter analogs and
around further away and older stars with current systems, and for future Giant
Segmented Mirror Telescopes to enable habitable exoplanet imaging[2].

In this paper we present three new AO and HCI technologies in ongoing
development in the Laboratory for AO (LAO) at UCSC, on the Santa Cruz Extreme
AO Laboratory (SEAL) testbed, which is dedicated to developing and testing
exoplanet imaging-focused HCI and ExAO technologies. SEAL in the context of
this paper is introduced in §2, the three SEAL projects are then discussed in
§3-5, and we then conclude in §6. A common theme to all three projects that we
will refer to throughout this paper is the use of “second stage AO” and/or
“multi-wavefront sensor (WFS) single conjugate AO (SCAO).” To clarify here,
“second stage AO” represents a “cascaded” AO system configuration, where a
second deformable mirror (DM) is placed downstream of the
beamsplitter/dichroic to the “first stage” WFS path (i.e., the second stage DM
is non-common path to the first stage WFS), with a separate second WFS driving
that second DM. Multi-WFS SCAO instead utilizes multiple WFSs pointing at the
same natural or laser guide star to control one or more common path DM(s),
which we will introduce and discuss further in §3.1. We will also note that
all laboratory results presented in this paper are planned for future
publications in peer reviewed journals, and so only a concise summary of
results at this time is presented in this paper.

##  2 SEAL LAYOUT

Figure 1 illustrates the optical layout of our SEAL testbed, previously
introduced in References 3 and 4.

Figure 1: Layout of the SEAL testbed, previously introduced in References 3
and 4. This figure (image credit: D. Dillon) is also additionally presented in
the following papers in this conference proceedings: References 5, 6, and 7.
SEAL currently operates at $\lambda=635$nm with HeNe lasers. There are
currently 4 DMs (SLM, segmented Iris AO DM, low order ALPAO DM, and high order
BMC kiloDM) and 6 WFS arms (SHWFS, PSF imaging camera, 3 sided Pyramid WFS,
Andor FAST WFS/4 sided PyWFS, Zernike WFS, vortex coronagraph arm). Most
hardware is controlled using a clone of the Keck II Pyramid WFS RTC[8], with
high level software written in Python.

Several other papers in this proceedings use SEAL: Ref. 9 for testing of
Linear Quadratic Gaussian control algorithms, Ref. 5 for testing a reflective
3-sided Pyramid WFS, Ref. 6 for testing a vector Zernike WFS to control an
IrisAO DM representing segmented telescope piston errors, and Ref. 7 for
development of a high order spatial light modulator (SLM) for simulating
laboratory atmospheric turbulence. In this paper, we will present three
additional projects using SEAL, introduced below along with their relevant
SEAL hardware components:

  * •

§3: Multi-WFS SCAO applied to FAST-SHWFS operations. This demonstration uses
the Fast Atmospheric Self-coherent camera (SCC) Technique (FAST)[10, 4] and
Shack Hartmann WFS (SHWFS) sensors to control our LO (LO) ALPAO 97 actuator DM
(hereafter referred to as the LODM, for low-order DM) and Boston Micromachines
Corporation (BMC) high order DM (34 actuators across the beam diameter,
hereafter referred to as the HODM, for high order DM, and in combination WT,
for “woofer-tweeter”).

  * •

§4: Focal plane wavefront sensing via a pupil chopping scheme, including first
an external optical chopper device (§4.1, hereafter referred to as version 1,
or v1) and second an internal DM-based approach (§4.2, hereafter v2), the
latter requiring no additional hardware beyond what is already present in a
conventional AO-fed imaging system. This demonstration uses the LODM,
controlled for v1 by a custom optical chopper and Andor camera and for v2 by
the FLIR Blackfly PSF imaging camera.

  * •

§5: Testing the bright pyramid WFS (bPWFS), a concept introduced in Ref. 11
and now tested further with SEAL in this paper. Here we use the LODM, custom
FPM, and Andor camera for bPWFS tests.

All of these above techniques present an interesting potential to further
develop with the above-described multi-WFS SCAO concept, which we discuss
further next in §3. We also note here that in this paper we utilize SEAL
calibration procedures described in Ref. 4, both specifically for FAST and
more generally for SEAL use (e.g., LODM And HODM voltage-to-wavefront
conversion).

##  3 FAST MULTI-WFS SCAO

###  3.1 Concept

Figure 2 illustrates the concept of multi-WFS SCAO, applied specifically to
our SEAL FAST-SHWFS WT system.

Figure 2: High-level illustration of how our SEAL SHWFS and FAST sensors,
through a series of different modal groups and temporal filters, will control
the LODM (11x11 ALPAO) and HODM (34x34 MEMS). See the text for more details.

This diagram builds on the concept initially presented in Ref. 10, but is now
adapted specifically to the SEAL testbed. As in Ref. 10, the main goal is to
enable steady-state closed-loop control with two on-axis WFSs controlling two
common-path DMs, which is done by modal groupings and temporal filters such
that the main diagram with all loops closed is consistent with the upper right
table in Fig. 2, giving AO control authority to only one spatio-temporal
regime for one WFS. An acquisition procedure is also illustrated in this
diagram, accounting for the context that the SHWFS must first close the loop
on atmospheric turbulence (step 1) and reach residual AO phases within the
FAST capture range for the FAST low and high order loops to then succesively
close (steps 2 and 3, respectively). Below we walk through the different real-
time steps shown in Fig. 2 between given WFS image readout and sending a
corresponding DM command.

  1. 1.

For each WFS we split modal groups between the LODM and HODM (and for the
three FAST LODM sub-groups, as we found in Ref. 4).

  2. 2.

Real-time WFS frames (SHWFS center of gravity, or CoG, values, FAST $I_{-}$
values) are converted into modal space via a Image-to-mode (I2M) matrix-vector
multiplication (MVM).

  3. 3.

A given SHWFS modal group can then be temporally filtered (there are no FAST
temporal filters here, instead giving FAST full control authority to all of
it’s accessible spatio-temporal regimes) on a real-time continuous buffer of
modal coefficients. We have identified equation 22 in Ref. 12 to enable such
high pass filtering.

  4. 4.

Modes for path $i$ are multiplied by a gain $g_{i}$. For modal gain
optimization (MGO), this can be a vector of different optimal gain values for
each mode.

  5. 5.

Modes are converted into DM commands via another MVM. Optionally adding data-
driven predictive wavefront control (PWFC), which we will discuss further
below in §3.3, adds another MVM. If no MGO and/or PWFC is used for FAST (i.e.,
which does not need temporal filtering), the whole line can be reduced to a
single real-time MVM.

  6. 6.

Lastly a simple type I integrator controller as shown enables closed-loop
operations on both DMs. More complex controllers could be implemented here
instead.

Also note that the speeds and filter cutoffs as shown in Fig. 2 are for only a
starting configuration; once we can demonstrate stable performance we can look
at changing WFS frame rates and corresponding temporal filter cutoffs to
simulate different guide star magnitudes, atmospheric conditions, etc.

###  3.2 SHWFS-WT Loop

Figure 3 shows telemetry results from real-time control between the SHWFS and
WT DMs.

Figure 3: 200 Hz SHWFS WT control of input full-scale atmospheric frozen-flow
turbulence. The LODM and HODM are both used to simulate input turbulence.

In this setup LODM modal bases include tip/tilt/focus (TTF) modes and
influence functions (IFs) with TTF removed (labeled “ALPAO IFs” in Fig. 3). We
found this separate TTF modal group was helpful in closed-loop to enable gain
and leak adjustments for a leaky integrator controller specific to only these
modes. HODM modal bases are Fourier modes between 5 and 10 cycles/pupil (c/p),
where the lower limit is set not to overlap with the ALPAO Nyquist limit and
the upper limit is set below the 23$\times$23 SHWFS lenslet array sampling.
Although this setup under-samples HODM-SHWFS Fried geometry, this 23$\times$23
subaperture modes is the only configuration available for this Thorlabs WFS-20
in high speed mode, which enables on-chip CoG processing to provide a stable
and high frame-rate readout at $\sim$200 Hz (with Thorlabs benchmarking
possible speeds up to 1 kHz), as showing in the left panel of Fig. 3.

We generate input atmospheric turbulence using the LODM and HODM via the WT
framework presented in Ref. 12, in particular that the spatial frequency
cutoff to switch from LODM to HODM turbulence is set well below the LODM
Nyquist limit. Note that this cutoff for injected DM-based turbulence is
different from the LODM control frequency cutoff, the latter instead
controlling out to the LODM Nyquist limit. We use an offline pre-calibrated
interaction matrix (IM) between LODM LO injected turbulence modes, HODM
injected turbulence modes, and SHWFS slopes in order to project a full power
law input atmospheric turbulence in HODM space onto LODM and HODM commands to
be applied in real-time via a single MVM.

As shown in Fig. 3, we are still adjusting and tuning this SHWFS-WT loop, as
AO residual phases are not yet at the Extreme AO (ExAO) levels we would like.
A 5 m/s frozen flow single-layer wind speed is used here, but with higher wind
speeds planned in the future with more faster loop speeds and a more optimized
control framework (e.g., implementing MGO and/or PWFC). However, our first
attempts at SHWFS-WT MGO showed that open and closed loop PSDs do not match
the theoretical rejection trasnfer function given the measured loop speed.
Since this measured loop speed is for our Python serial implementation (i.e.,
measured times include WFS image acquisition, computational lag, and applying
DM commands), it is hard to identify what factors might be causing limited
performance gains. We have already started developing a low-level C-based
parallelized loop code that should help improve stability and performance
before integrating this loop with the FAST WT loop.

###  3.3 FAST WT Loop

Figure 4 shows closed-loop telemetry results for FAST running real-time
correction on HODM-injected AO residual turbulence.

Figure 4: 100 Hz FAST WT control in input frozen flow AO-residual atmospheric
turbulence. The HODM is used to simulate input turbulence. A 50 Hz loop speed
version of this figure was initially published in Ref. 4.

The panels are similar to Fig. 3, but with slightly different modal groups and
acquisition procedures. A 1 m/s frozen flow wind speed is used, with planned
higher wind speeds accessible by higher frame rates and/or more advanced real-
time control algorithms. These results are consistent with the FAST findings
in Ref. 4, where we demonstrated how four different modal groups enabled good
linearity and closed-loop performance: TTF, LO Zernike modes (labeled
“Zernikes”), LO Fourier modes (labeled “LO Fourier,” still diffracting sine
spots within the FAST focal plane mask inner working angle), and high order
Fourier modes (labeled “HO Fourier,” using the conventional modal basis to
generate a half dark hole, or DH). The former three modes are controlled with
the LODM, where as the latter high order Fourier modal group is controlled
with the HODM. We also implement the same LO acquisition procedure as in Ref.
4 and as in Fig. 2: first closing TTF, then closing LODM Zernike and Fourier
modes, and lastly closing HODM Fourier modes, which we found was necessary to
enable stable closed-loop DHs over time. Again, similar to §3.2, we are still
tuning parameters for optimal closed-loop, including ongoing developments to
migrate our Python real-time control software into C in a format consistent
with Fig. 2.

We have also demonstrated real-time FAST PWFC, to our knowledge a first of
it’s kind for focal plane WFSs, shown in Fig. 5.

(a) (b)

Figure 5: FAST PWFC for LO Zernike modes of frozen flow AO residual
turbulence, showing the time domain version in (a), running at $\sim$50 Hz and
the time-collapsed histograms in (b). The shaded regions in each panel
represent the standard deviation around the mean represented by the line of
the same color.

Here we are taking a similar approach to what our UCSC LAO group is currently
testing at W. M. Keck Observatory[13], we implement PWFC using Empirical
Orthogonal Functions (EOF)[14] method. Using EOF we predict the closed-loop
low order modal coefficients (note that this is different from the original
EOF implementation that requires a pseudo open loop reconstruction) by storing
up modal coefficients to generate a predictive filter. The filter is then used
in a MVM to subsequently predict the FAST coefficients over the system delay,
represented by the “M2A” boxes in Fig. 2. In Fig. 5, we show preliminary
results using EOF with FAST for LO modes. The distributions of residual
wavefront error (WFE) over time indicate an improvement of 1.3x when EOF
(i.e., PWFC) is used, with a WFE standard deviation improvement of 1.9x (the
latter indicating increased stability). We plan to further improve the
performance of FAST PWFC by optimizing the predictor order and the amount of
training data.

##  4 PUPIL CHOPPING FOR FOCAL PLANE WAVEFRONT SENSING

We present two versions of the pupil chopping concept here: the first version
(hereafter v1) in §4.1 uses an external optical chopping device. The second
version (hereafter v2) in §4.2 is a newer concept that continues to be
developed in simulation and tested in the lab (project lead: co-I J. Perez-
Soto), and so we present preliminary results in §4.2.

###  4.1 v1: External Optical Chopping

Figure 6 illustrates the concept of focal plane wavefront sensing with a pupil
plane optical chopper.

(a) (b)

Figure 6: Illustration of the pupil chopping concept, both from a geometric
(a) and Fourier optics (b) perspective. In (b), all focal plane images are
shown with the same field of view. Panels 2-3 are shown on a log scale; all
other panels are shown on a linear scale. Zernike and Fourier modes in panels
5-12 all use 1 nm amplitudes. Panels 7-12 indicate that pupil chopping for
focal plane wavefront sensing can resolve the sign ambiguity of focus and the
relative phase of Fourier modes, which is not possible with a single PSF
image.

In short, two focal plane images are required to enable a WFS measurement: one
with the chopper blade partially blocking the pupil, and one with the pupil
unblocked; the latter can be used for science while the former is limited to
wavefront sensing. A “reference” WFS frame is saved before on-sky operations
and then used analogously to the concept of on-sky reference slopes for other
pupil plane WFSs, enabling a real-time WFS frame to in principle flatten the
wavefront to the level of the pre-calibrated reference WFS frame. It should be
noted that this technique is complementary to pupil amplitude diversity-based
absolute phase retrieval methods, such as the asymmetric pupil Fourier
WFS[15]. This technique enables a linear-least squares reconstruction, optimal
for high speed wavefront control of AO residuals but not able to flatten the
absolute phase below a pre-calibrated best flat (e.g., due to evolving on-sky
quasi-static aberrations), whereas phase retrieval methods can use the chopped
pupil PSF image to recover the absolute phase to track such quasi-static
effects, but typically requiring non-linear iterative algorithms that cannot
be run at high speed.

Our SEAL laboratory setup for v1 pupil chopping presented in this subsection
uses the LODM and Andor camera as previously descrived, and additionally a
Thorlabs MC2000B optical chopper with blade MC1F10, and a digital signal
processing controller and custom electrical doubler so that the Andor camera
readout is synchronized with the optical chopper (with the chopper as the
leader and the Andor camera as the follower111Although the AO field has
historically adopted “master and slave” terminology from the computer science
field to refer to the chopper and camera in this configuration, respectively,
it is now clear that such discriminatory terminology should no longer be used.
Following the computer science field, instead we will refer to such
configuration as “leader and follower” throughout this paper, encouraging
other members of the AO field to do the same.) to produce a chopped and un-
chopped pupil every other frame, respectively.

Using a Zernike modal basis of up to five radial orders (limited to higher
orders by phase jitter of the chopper blade in our lab setup), on-air closed-
loop results are shown in Fig. 7.

(a) (b)

Figure 7: Real-time correction of on-air wavefront error (at a frame rate of
about 25 Hz), showing both tip/tilt/focus (TTF) and all controlled modes
($n=1\rightarrow 5$), both in the time domain (a) and by WFE distribution
comparing the open and closed loop sequences (b; i.e., respectively to the
left and right of the dotted line in panel a). Note that the y axis in panel
(a) and x axis in panel (b) are not correctly normalized; only the relative
reduction between open and closed loop WFE should be considered for the
purposes of this paper.

Although the highly stabilized environment from our granite optical table
allows only a moderate gain, Fig. 7 clearly demonstrates a closed-loop WFE
reduction relative to open loop for all controlled modes with only on-air
temporal disturbances present, illustrating the potential of this technique to
stabilize quasi-static temporal disturbances on-sky (e.g., thermal variations,
flexure, and/or other sources of slow beam wander). Note these results help
address concerns that the spinning chopper wheel would be generating
additional on-air turbulence, since telemetry from closing the loop on-air
clearly shows improvement compared to open-loop.

Figure 8: Open and closed loop image stacks of un-chopped images (left and
middle panels) and corresponding radial intensity profile (right panel) for
frozen flow 10 m/s AO residual turbulence applied on the DM in real time at
$\sim$50 Hz. The right panel shows a $\sim$20% Strehl improvement by closing
the loop on AO residuals.

Building on these on-air results, Fig. 8 shows un-chopped PSF images and a
corresponding radial profile for LODM-injected AO residual turbulence with and
without closing the loop on that turbulence, clearly illustrating the
potential of this technique for second stage AO and/or multi-WFS SCAO. In
addition to the above-mentioned chopper blade phase jitter limiting our
ability to measure and control Zernike modes beyond 5 radial orders (for which
a custom chopper solution could reduce sufficiently to enable higher order
control, such as Ref. 16), a more fundamental complexity is that an additional
optical chopper hardware component must be designed and installed into
instruments for this technique to work, increasing project cost and
complexity. Science duty cycle with this v1 technique is also limited to, at
most 50%. Instead, v2, discussed next in §4.2, improves on all of these issues
(i.e., removes the phase jitter problem, does not require additional optical
chopper hardware, and increases the science duty cycle).

###  4.2 v2: DM-based Chopping

Instead of using a physical optical chopper device, v2 of this approach uses
the DM of an existing AO system to apply a local tilt to a fraction of the
illuminated pupil. PSF images are recorded where this DM “chop” commands are
added to all other DM commands every other frame. Fig. 9 illustrates this
concept using laboratory images from SEAL.

Figure 9: DM-based pupil chopping concept illustration, using SEAL testbed
images, clearly demonstrating the potential of this to resolve the sign of
focus and other common degenerates present from a single PSF image, but still
enable high duty cycle science imaging.

At the moment this local tilt uses a peak-to-valley stroke of $\sim$1$\mu$m
with our LODM, but we are still investigating whether lower stroke options
(e.g., enabling HODM-based pupil chopping without saturating available DM
stroke) are compatible with this technique. Although Fig. 9 only shows low
order wavefront modes, this technique is broadly applicable to wavefront
sensing and control of spatial orders out to the DM Nyquist limit, and by
being a focal plane wavefront sensor it benefits from natural wavefront
spatial filtering properties, enabling binary masks on the focal plane image
to act as a natural anti-aliasing filter and minimize non-linear cross talk
between modes (e.g., as demonstrated for FAST in Ref. 4). Using the focal
plane image and such binary masks for wavefront sensing with this technique
also distinguishes this technique from the differential optical transfer
function method[17].

Figure 10 shows preliminary DM-based pupil chopping linearity simulations for
low order modes.

Figure 10: LO linearity simulations for our proposed DM-based pupil chopping
technique, assuming $\lambda$=1.6$\mu$m to show peak-to-valley (PV) units on
both axes.

As expected from Fig. 9, LO modes show excellent linearity for phases
$\lesssim$1 rad rms, showing applicability to correct AO residual phases but
not open loop atmospheric turbulence. Further simulations and SEAL experiments
with this DM-based pupil chopping technique are ongoing.

The hardware setup of this technique is extremely simple, requiring only a PSF
imaging camera in addition to existing AO system hardware, making this
compatible with most current and future AO systems. A high-speed PSF imaging
camera would enable second stage and/or multi-WFS SCAO correction of AO
residuals, but even a slower readout camera could enable on-sky stabilization
of quasi-static wavefront errors. This technique also has broad applicability
to other AO fields (e.g., laser guide stars, off-axis stars, crowded star
fields, extended sources); as long as a camera images the focal plane of the
intended light source for guiding and the DM is conjugated to the pupil plane,
the principles of this technique remain unchanged. Increasing the science duty
cycle above 50% is also feasible with this setup (e.g., if the average PSF
camera readout is at 1 kHz, if the guide star is bright enough the chopped
pupil frame can be a 0.1 ms exposure and the un-chopped pupil science frame
can be 0.9 ms, enabling a 90% science duty cycle). However, a multi-WFS SCAO
approach with this technique means that the other AO WFS sees the local DM
tilt every other frame. Possible solutions to this are to (1) remove this
local tilt mode from the control basis for the first stage WFS, (2) drop every
frame when the DM has this local tilt applied (which if operating, e.g., at a
90% duty cycle should yield a minimal performance loss), or (3) use this
technique only for second stage AO, meaning that only a second DM that is not
common path to the first stage WFS is driven by this DM-based pupil chopping
technique. All options should be investigated further to more quantitatively
understand the benefits and disadvantages of each approach.

##  5 TESTING THE BRIGHT PYRMAID WFS

Lastly, we have recently been using SEAL to further test and develop the
bright Pyramid WFS (bPWFS) concept, initially presented in Ref. 11 and
illustrated here in Fig. 11.

(a) (b)

Figure 11: Illustration of the bPWFS concept, first introduced in Ref. 11. In
panel a, we show fabricated masks (i.e., Zygo interferometer-measured depth
profiles) from KAUST and corresponding pupil images from SEAL after aligning
these masks, illustrating that the non-modulated bPWFS diffracts more light
into the pupil footprints than the non-modulated PWFS (the PWFS and bPWFS
pupil images are both normalized to the same scale). The benefits of that
light redistribution are illustrated in panel b, showing simulations (adapted
from Ref. 11 to integrate over spatial frequencies less than 16 c/p for an
assumed 32 actuator diameter DM), showing that for input AO residual phases,
the bPWFS vs. PWFS output improves the non-linearity, aliasing, and photon
noise error budget components.

In this paper, we present further testing of the bPWFS concept on SEAL, using
custom-made reflective 16-level scalar Aluminum masks fabricated by the
Computational Imaging group at King Abdullah University of Science and
Technology (KAUST). This fabrication process, reaching $\sim$3 micron spatial
resolution, $\sim$few nm-level depth precision, but limited to $\sim$800 nm
dynamic range (hence the need to phase wrap the PWFS tilt angles), is
described in detail in Ref. 18. Despite the need for phase wrapping, which is
inherently monochromatic, the high depth precision and spatial resolution
could potentially outweigh chromaticity effects from phase wrapping to improve
PWFS quality compared to existing fabrication techniques, and would be
interesting to consider for applications in other PWFS systems. In Fig. 11,
our high quality KAUST PWFS and bPWFS masks made specifically for SEAL clearly
show more light within bPWFS vs. PWFS pupil footprints, illustrating the
potential to demonstrate enhanced performance in the lab as simulations
suggest.

Our first approach was to consider a closed-loop bootstrapping procedure where
a modulated bPWFS is first used to close the loop on atmospheric turbulence
and then modulation is switched off to see the non-modulated bPWFS gains.
However, we found in simulations, both for the PWFS and bPWFS, that non-
linearities dominate the error budget when modulation is switched off and as
such performance immediately degrades. Instead, in this paper we only
investigate non-modulated linearities for the two WFSs, keeping in mind that
future solutions will be needed to integrate this non-modulated bPWFS
performance gain into real AO systems, such as a second stage bPWFS, multi-WFS
SCAO approach, and/or non-linear bPWFS reconstruction algorithms.

Figure 12 shows both simulations, using OOMAO[19], and corroborating SEAL
testbed results that the non-modulated bPWFS is indeed more linear than the
non-modulated PWFS, both for low and high order Zernike modes.

Figure 12: OOMAO simulations and corroborating SEAL data (left and right
columns, respectively) showing the non-modulated bPWFS is more linear (i.e.,
both linear and capture range) than the regular PWFS for both low and high
order Zernike modes (top and bottom rows, respectively). Due to alignment
issues, SEAL LO linearity compares the integral of the absolute value of
linearity curves for different LO Zernike modes, showing the ratio of this
integral for the bPWFS relative to the PWFS (i.e., a number greater than one
indicates enhanced bPWFS linearity). This cumulative linearity metric works
for our lab setup where there is only small misalignment between the PWFS and
bPWFS masks, but it would not be representative for larger alignment errors
where the probed linear range does not reach the saturated region on both
positive and negative ends of a given linearity curve.

Interestingly, one astigmatism mode shows no linearity gain in simulation and
a slight degradation in lab data, which should be further investigated. The
choice of modal basis also needs further investigation; initially we used TTF
and then LODM actuator pokes with TTF removed (i.e., as done for LODM SHWFS
control in §3.2), but simulations and lab data showed negligible bPWFS
linearity enhancement with this basis. Particularly for lab data, we noticed
that PWFS IM images may have had insufficient signal-to-noise ratios for the
small actuator poke amplitudes needed to be within the linear regime, but
where on-air turbulence averaging over longer sequences of stacked images also
degraded IM quality; the optimal modal basis that can be used to acquire IMs
at high enough speed to freeze on-air turbulence should there be further
investigated. Alignment also served as a potential limitation in laboratory
demonstrations. Although we wrote an automated alignment script designed to
minimize the standard deviation of the pupil intensities over tip, tilt, and
focus to find the best alignment in an automated way each time we swapped PWFS
and bPWFS masks on SEAL, LO linearity curves remained offset between the two
masks (hence the need to show their integrals in Fig. 12 for a fair
comparison); more robust alignment algorithms (e.g., alignment by slope
minimization) could enable a more uniform comparison between the two masks.

As discussed above, future work will need to include a more integrated plan on
how to operate the non-modulated bPWFS in closed-loop within existing AO
systems. Before that, future bPWFS SEAL testing goals not already mentioned
include linearity analysis of higher order modes beyond the LODM Nyquist limit
with our HODM and closed-loop analysis of AO residual level input WFE. For the
latter, our SLM[7] can be used for an aliasing error budget analysis,
comparing DM-injected and reconstructed turbulence can be used for a non-
linearities error budget analysis, and results with our tuneable light source
(with a range from 0.01 to 8 mW) comparing a photon noise-limited and WFE-
limited regime can be used for a photon noise error budget analysis.

##  6 CONCLUSIONS AND FUTURE WORK

In this paper we have presented laboratory results with SEAL (§2) on three
ongoing areas of HCI and ExAO technology development:

  1. 1.

§3: Multi-WFS SCAO using a SHWFS and FAST to control a LODM and HODM. Thus far
we have demonstrated separate high-speed real-time control of each WFS with
both DMs but have not yet combined both WFSs into a single real-time loop.
SEAL development to implement real-time temporal filtering in order to combine
both WFSs is ongoing.

  2. 2.

§4: Pupil chopping, both with an external optical chopper (§4.1) and internal
DM-based chopping approach (§4.2). We have demonstrated the power of either
approach to enable second stage AO and/or multi-WFS SCAO, in this case for a
non-coronagraphic science image at a $>$50% science duty cycle (but with
coronagraphic pupil chopping investigations to come in the future).

  3. 3.

§5: testing the bright Pyramid WFS concept. We have demonstrated that the non-
modulated bPWFS is more linear than the non-modulated PWFS, but also finding
that a modulated to non-modulated closed-loop bootstrapping procedure is
ineffective, prompting further work needed on the feasibility of second stage
and/or multi-WFS SCAO with a bPWFS.

###### Acknowledgements.

We gratefully acknowledge research support of the University of California
Observatories and UCSC’s Lamat NSF REU program for funding this research. This
work also benefited from the 2022 Exoplanet Summer Program in the Other Worlds
Laboratory (OWL) at the University of California, Santa Cruz, a program funded
by the Heising-Simons Foundation. Author B. Gerard thanks the 2018 SCExAO team
for hosting discussions that led to the pupil chopping concept and the UCSC
LAO group, particularly Dominic Sanchez, for discussions that led to the DM-
based pupil chopping idea.

## References

  * [1] Nielsen, E. L., De Rosa, R. J., Macintosh, B., Wang, J. J., Ruffio, J.-B., Chiang, E., Marley, M. S., Saumon, D., Savransky, D., Ammons, S. M., Bailey, V. P., Barman, T., Blain, C., Bulger, J., Burrows, A., Chilcote, J., Cotten, T., Czekala, I., Doyon, R., Duchêne, G., Esposito, T. M., Fabrycky, D., Fitzgerald, M. P., Follette, K. B., Fortney, J. J., Gerard, B. L., Goodsell, S. J., Graham, J. R., Greenbaum, A. Z., Hibon, P., Hinkley, S., Hirsch, L. A., Hom, J., Hung, L.-W., Dawson, R. I., Ingraham, P., Kalas, P., Konopacky, Q., Larkin, J. E., Lee, E. J., Lin, J. W., Maire, J., Marchis, F., Marois, C., Metchev, S., Millar-Blanchaer, M. A., Morzinski, K. M., Oppenheimer, R., Palmer, D., Patience, J., Perrin, M., Poyneer, L., Pueyo, L., Rafikov, R. R., Rajan, A., Rameau, J., Rantakyrö, F. T., Ren, B., Schneider, A. C., Sivaramakrishnan, A., Song, I., Soummer, R., Tallis, M., Thomas, S., Ward-Duong, K., and Wolff, S., “The Gemini Planet Imager Exoplanet Survey: Giant Planet and Brown Dwarf Demographics from 10 to 100 au,” Astronomical Journal 158, 13 (July 2019).
  * [2] Jensen-Clem, R., Hinz, P., Skemer, A., Wizinowich, P., Jovanovic, N., Mazin, B. A., John I. Bailey, I., Frazin, R. A., Sallum, S., Males, J. R., and Tamura, M., “A Technology and Science Gap List for Habitable-Zone Exoplanet Imaging with Ground-Based Extremely Large Telescopes,” in [In this proceeding ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (July 2022).
  * [3] Jensen-Clem, R., Dillon, D., Gerard, B., van Kooten, M. A. M., Fowler, J., Kupke, R., Cetre, S., Sanchez, D., Hinz, P., Laguna, C., Doelman, D., and Snik, F., “The Santa Cruz Extreme AO Lab (SEAL): design and first light,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11823, 118231D (Sept. 2021).
  * [4] Gerard, B. L., Dillon, D., Cetre, S., and Jensen-Clem, R. M., “Laboratory demonstration of real-time focal plane wavefront control of residual atmospheric speckles,” Journal of Astronomical Telescopes, Instruments, and Systems 8(3), 1 – 22 (2022).
  * [5] Sanchez, D. F., Hinz, P., and Dillon, D., “Preliminary lab demonstration of a 3-sided reflective pyramid wavefront sensor for Shane AO using SEAL testbed,” in [In this proceeding ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (July 2022).
  * [6] Salama, M., Jensen-Clem, R., van Kooten, M., Dillon, D., Gerard, B. L., Fowler, J., Cetre, S., Snik, F., and Doelman, D., “Vector Zernike Wavefront Sensor on the Santa Cruz Extreme AO Lab (SEAL) Testbed,” in [In this proceeding ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (July 2022).
  * [7] van Kooten, M. A., Jensen-Clem, R., Fowler, J., Dillon, D., Kupke, R., Salama, M., and Gerard, B. L., “Spatial light modulator on Santa Cruz Extreme AO Laboratory (SEAL) testbed,” in [In this proceeding ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (July 2022).
  * [8] Cetre, S., Guyon, O., Bond, C., Chun, M., Mawet, D., Wizinowich, P., Lockhart, C., Goebel, S., and Wetherell, E., “A near-infrared pyramid wavefront sensor for Keck adaptive optics: real-time controller,” in [Adaptive Optics Systems VI ], Close, L. M., Schreiber, L., and Schmidt, D., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10703, 1070339 (July 2018).
  * [9] Sengupta, A. R., Gerard, B. L., Dillon, D., van Kooten, M., Gavel, D., and Jensen-Clem, R., “Laboratory Demonstration of Optimal Identification and Control of Tip-Tilt Systems,” in [In this proceeding ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (July 2022).
  * [10] Gerard, B. L., Véran, J.-P., Singh, G., Herriot, G., Lardière, O., and Marois, C., “Fast focal plane wavefront sensing as a second stage adaptive optics wavefront sensor,” in [Adaptive Optics Systems VII ], Schreiber, L., Schmidt, D., and Vernet, E., eds., 11448, 483 – 499, International Society for Optics and Photonics, SPIE (2021).
  * [11] Gerard, B. L., Chambouleyron, V., Jensen-Clem, R., and Sauvage, J.-F., “The bright pyramid wavefront sensor,” in [Techniques and Instrumentation for Detection of Exoplanets X ], Shaklan, S. B. and Ruane, G. J., eds., 11823, 403 – 412, International Society for Optics and Photonics, SPIE (2021).
  * [12] Gavel, D. and Norton, A., “Woofer-tweeter deformable mirror control for closed-loop adaptive optics: theory and practice,” in [Adaptive Optics Systems IV ], Marchetti, E., Close, L. M., and Vran, J.-P., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9148, 91484J (Aug. 2014).
  * [13] van Kooten, M. A. M., Jensen-Clem, R., Cetre, S., Ragland, S., Bond, C. Z., Fowler, J., and Wizinowich, P., “Predictive wavefront control on Keck II adaptive optics bench: on-sky coronagraphic results,” Journal of Astronomical Telescopes, Instruments, and Systems 8, 029006 (Apr. 2022).
  * [14] Guyon, O. and Males, J., “Adaptive Optics Predictive Control with Empirical Orthogonal Functions (EOFs),” arXiv e-prints , arXiv:1707.00570 (July 2017).
  * [15] Martinache, F., “The Asymmetric Pupil Fourier Wavefront Sensor,” Publications of the Astronomical Society of the Pacific 125, 422 (Apr. 2013).
  * [16] Johnson, A. B., Marois, C., Gamroth, D., Lardière, O., Thompson, W., Singh, G., Fitzsimmons, J., and Bradley, C., “Blinking the fringes, initial development and results of the Ultra-Low Speed Optical Chopper for the Self-Coherent Camera,” in [In this proceeding ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (July 2022).
  * [17] Codona, J. L., “Differential OTF Wavefront Sensing,” arXiv e-prints , arXiv:1309.6719 (Sept. 2013).
  * [18] Fu, Q., Amata, H., Gerard, B., Marois, C., and Heidrich, W., “Additive lithographic fabrication of a Tilt-Gaussian-Vortex mask for focal plane wavefront sensing,” in [Optifab 2021 ], Nelson, J. D. and Unger, B. L., eds., 11889, 162 – 170, International Society for Optics and Photonics, SPIE (2021).
  * [19] Conan, R. and Correia, C., “Object-oriented Matlab adaptive optics toolbox,” in [Adaptive Optics Systems IV ], Marchetti, E., Close, L. M., and Vran, J.-P., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9148, 91486C (Aug. 2014).


# Inverse properties of a class of pentadiagonal matrices related to higher
order difference operators

Bakytzhan Kurmanbek  Nazarbayev University, Department of Mathematics, 53
Kabanbay Batyr Ave, Nur-Sultan 010000, Kazakhstan  Yogi Erlangga  Zayed
University, Department of Mathematics, Abu Dhabi Campus, P.O. Box 144534,
United Arab Emirates  Yerlan Amanbek  Nazarbayev University, Department of
Mathematics, 53 Kabanbay Batyr Ave, Nur-Sultan 010000, Kazakhstan
Corresponding author

###### Abstract

This paper analyzes the convergence of fixed-point iterations of the form
$\bm{u}=f(\bm{u})$ and the properties of the inverse of the related
pentadiagonal matrices, associated with the fourth-order nonlinear beam
equation. This nonlinear problem is discretized using the finite difference
method with the clamped-free and clamped-clamped boundary conditions in the
one dimension. Explicit formulas for the inverse of the matrices and norms of
the inverse are derived. In iterative process, the direct computation of
inverse matrix allows to achieve an efficiency. Numerical results were
provided.

Keywords. explicit formula, pentadiagonal matrices, finite difference,
nonlinear beam equation, fixed point method

##  1 Introduction

Many applications give arise to mathematical problems that involve numerical
computations with pentadiagonal matrices, which require their inversion (see
[1] and references therein). Even though inversion of a nonsingular
pentadiagonal matrix can be done efficiently by a numerical linear algebra
software, explicit inverse formulas are useful, for example, in a computer
algebra software.

Early results on inverses of banded matrices can be traced as far back as to
the work of [2, 3, 4] for general band matrices. Results for band Toeplitz
matrices are given in [5], with explicit inverse formulas for tridiagonal
matrices in [6, 7] and pentadiagonal matrices in [8, 9, 10, 1, 11, 12, 13]. In
addition, properties including determinants of such matrices related to finite
difference operators have been investigated, e.g. in [14, 15, 16].

In this study, we focus on the specific pentadiagonal matrices arising in a
fixed-point iteration for numerically solving the fourth-order nonlinear beam
equation:

$\displaystyle\displaystyle\frac{d^{4}\widehat{\phi}}{d\widehat{x}^{4}}=\alpha_{1}e^{\displaystyle-\alpha_{2}\widehat{\phi}},\quad\widehat{x}\in\Omega=(0,L).$

This nonlinear equation finds applications in mechanical and civil
engineering, which models, e.g., a cantilever beam subjected to swelling
pressure on one side. In the above equation, the right-hand side term is the
swelling pressure, which in this form is proposed by Grob [17], based on
empirical studies (see, e.g., [18] and the references therein), $L>0$ is the
length of the beam, and $\alpha_{1},\alpha_{2}>0$ represents the mechanical
property of the beam, which are assumed to be constant.

Scaling the domain to unity using the dimensionless variable $x=\widehat{x}/L$
and setting $\phi=\alpha_{2}\widehat{\phi}$ yields

$\displaystyle\displaystyle\frac{d^{4}\phi}{dx^{4}}=Ke^{-\phi},\text{ in
}\Omega=(0,1),$ (1)

where $K=\alpha_{1}\alpha_{2}L>0$. We shall use this formulation throughout.
For (1) two types of boundary conditions are employed:

  1. 1.

Clamped-Free (CF) condition:

$\displaystyle\phi(0)=\phi^{\prime}(0)=0\quad\text{and}\quad\phi^{\prime\prime}(1)=\phi^{\prime\prime\prime}(1)=0,$
(2)

  2. 2.

Clamped-Clamped (CC) condition:

$\displaystyle\phi(0)=\phi^{\prime}(0)=0\quad\text{and}\quad\phi(1)=\phi^{\prime}(1)=0.$
(3)

Since $\displaystyle\frac{d^{4}\phi}{dx^{4}}=Ke^{-\phi}>0$, obviously,
$\phi=0$ can not be a solution, even though it satisfies the boundary
conditions.

The solution of (1) with the boundary conditions (2) is concave up and an
increasing function, which can be deduced from a mixed formulation of (2):

$\begin{cases}\displaystyle\frac{d^{2}\omega}{dx^{2}}=Ke^{-\phi},&\omega(1)=\omega^{\prime}(1)=0,\\\\[10.0pt]
\displaystyle\frac{d^{2}\phi}{dx^{2}}=\omega,&\phi(0)=\phi^{\prime}(0)=0.\\\
\end{cases}$ (4)

From the first part of (4), with $e^{-\phi}>0$ in $\Omega$,
$\omega^{\prime\prime}>0$, and $w^{\prime}$ increases in $\Omega$. The
condition $\omega^{\prime}(1)=0$ requires that $w^{\prime}<0$ in $\Omega$,
which furthermore, together with the condition $\omega(1)=0$, implies that
$\omega>0$ and decreases. From the second part of (4), we have
$\phi^{\prime\prime}=\omega>0$; thus, $\phi^{\prime}$ is an increasing
function in $\Omega$. Since $\phi^{\prime}(0)=0$, $\phi>0$, which implies
$\phi>0$ and increases. This characterization also holds in the finite-
difference setting based on the second-order scheme we use in this paper (c.f.
Section 4).

Numerical methods based on finite element methods for (1) have been proposed
and studied, e.g., in [19, 20], where focus is given on the accurate
approximation of the solution. This paper approaches the problem from a
different angle, with emphasis put on the convergence of the iteration method
of the form

$\phi=\mathcal{L}^{-1}\left(Ke^{-\phi}\right),$

where $\mathcal{L}=d^{4}/dx^{4}$, and the properties of the related iteration
matrices involved. Using the second-order finite difference approach, these
matrices are pentadiagonal and near Toeplitz.

In this paper, we present explicit formulas for inverses of the specific
pentadiagonal matrices and their bounds of norms, which are necessary in the
convergence analysis of the fixed-point iteration. As the inverse can be
formed explicitly, we are able to construct an exact norm of some of those
matrices. The convergence rate for the clamped-free and clamped-clamped
problems were derived and then numerical examples were presented for different
parameters.

The paper is organized as follows. Section 2 is devoted to the convergence and
the inverse of the iteration matrix for problem with the clamp-free condition.
Similar discussion for the clamp-clamp condition is given in Section 3.
Numerical results are presented in Section 4, followed by some concluding
remarks in Section 5.

##  2 The case with clamped-free boundary conditions

We consider $n+1$ equidistant grid points on the closed interval $[0,1]$, with
the distance (grid size) $h=1/n$, at which the solution of (1) is approximated
by a finite difference scheme. Each grid point is indexed by $i=0,\dots,n$,
where $i=0$ and $n$ correspond to the boundary points. Throughout the paper,
we shall consider $n\geq 5$ for $A$ to be a meaningful approximation to the
differential operator $\mathcal{L}$, even though $n=5$ may not be of practical
interest.

For the interior nodes, $1\leq i\leq n-1$, the fourth-order derivative is
approximated by the second-order finite difference scheme:

$\frac{d^{4}\phi}{dx^{4}}(x_{i})\approx\frac{1}{h^{4}}(\phi_{i-2}-4\phi_{i-1}+6\phi_{i}-4\phi_{i+1}+\phi_{i+2}),$

where $x_{i}=ih$ and $\phi_{i}\equiv\phi(x_{i})$. For $i=2$, we just impose
the boundary condition $\phi(0)\equiv\phi_{0}=0$. For $i=1$, $\phi_{-1}$
corresponds to a fictitious point outside the computational domain, which is
eliminated using the central scheme approximation to the boundary condition
$\phi^{\prime}(0)=0$. Similar approaches are used for $i=n-1$ and $n$, with
the boundary conditions
$\phi^{\prime\prime}(1)=\phi^{\prime\prime\prime}(1)=0$ be approximated by
appropriate second-order finite difference schemes.

The resultant system of nonlinear equations is

$\displaystyle A\bm{u}=h^{4}K\exp(\bm{-u}),$ (5)

where $\bm{u}=(u_{1},\dots,u_{n})^{T}\in\mathbb{R}^{n}$, with
$u_{i}\approx\phi(x_{i})$, and

$\displaystyle A:=\left[\begin{array}[]{rrrrrrr}7&-4&1&0&\cdots&&0\\\
-4&6&-4&\ddots&\ddots&&\vdots\\\ 1&-4&\ddots&\ddots&&&\\\
0&\ddots&\ddots&&\ddots&1&0\\\ &\ddots&&\ddots&6&-4&1\\\ \vdots&&&1&-4&5&-2\\\
0&\cdots&&0&2&-4&2\end{array}\right].$ (13)

Here, $A\in\mathbb{R}^{n\times n}$ is a nonsymmetric, nondiagonally dominant
pentadiagonal matrix.

Our first result on $A$ is that it is nonsingular. In fact, we have the
following theorem of the explicit inverse of matrix

######  Theorem 2.1.

Let $B=[b_{i,j}]_{i,j=1,n}\in\mathbb{R}^{n\times n}$ such that

$\displaystyle b_{i,j}$
$\displaystyle=\displaystyle\frac{3ij^{2}+j-j^{3}}{6},\quad\forall j\leq
i,i\in\\{1,2,\dots,n\\},j\in\\{1,2,\dots,n-1\\},$ $\displaystyle b_{i,n}$
$\displaystyle=\frac{1}{2}b_{n,i},$ $\displaystyle b_{n,n}$
$\displaystyle=\frac{1}{12}n(2n^{2}+1),$ $\displaystyle b_{i,j}$
$\displaystyle=b_{j,i},\quad i,j\in\\{1,2,\dots,n-1\\}.$

Then $B$ is the inverse of $A$, where $A=[a_{i,j}]_{i,j=1,n}$ is given in
(13).

###### Proof.

The proof is done by the direct computation. Let $D$ be matrix such that
$D=AB$. We want to show that the product

$d_{i,j}:=\left[a_{i,1}\,\,a_{i,2}\,\,\cdots\,\,a_{i,n}\right]\left[\begin{array}[]{c}b_{1,j}\\\
b_{2,j}\\\ \vdots\\\ b_{n,j}\end{array}\right]=\begin{cases}1,&i=j,\\\
0,&i\neq j.\end{cases}$

In other words, $D$ is the identity matrix $n\times n$.

  1. (i)

The case $3\leq i\leq n-2$ and $1\leq j\leq n$.

In this case, $a_{i,i-2}=1$, $a_{i,i-1}=-4$, $a_{i,i}=6$, $a_{i,i+1}=-4$,
$a_{i,i+2}=1$, while the others are $0$. Therefore,

$d_{i,j}=b_{i-2,j}-4b_{i-1,j}+6b_{i,j}-4b_{i+1,j}+b_{i+2,j}.$ (14)

If $i=j$, then $b_{i-2,i}=b_{i,i-2}=(2i^{3}-6i^{2}+i+6)/6$,
$b_{i-1,i}=b_{i,i-1}=(2i^{3}-3i^{2}+i)/6$, $b_{ii}=(2i^{3}+i)/6$,
$b_{i+1,i}=(2i^{3}+3i^{2}+i)/6$, and $b_{i+2,i}=(2i^{3}+6i^{2}+i)/6$, yielding
$d_{i,i}=1$.

For $i\neq j$, we consider several cases.

    1. (a)

$j\leq i-2$; Then $b_{i-2,j}=(3ij^{2}+j-6j^{2}-j^{3})/6$,
$b_{i-1,j}=(3ij^{2}+j-3j^{2}-j^{3})/6$, $b_{i,j}=(3ij^{2}+j-j^{3})/6$,
$b_{i+1,j}=(3ij^{2}+j+3j^{2}-j^{3})/6$,
$b_{i+2,j}=(3ij^{2}+j+6j^{2}-j^{3})/6$, yielding $d_{i,j}=0$.

    2. (b)

$j=i-1$; Then $b_{i-2,j}=(2i^{3}-9i^{2}+13i-6)/6$,
$b_{i-1,j}=(2i^{3}-6i^{2}+7i-3)/6$, $b_{i,j}=(2i^{3}-3i^{2}+i)/6$,
$b_{i+1,j}=(2i^{3}-5i+3)/6$, $b_{i+2,j}=(2i^{3}+3i^{2}-11i+6)/6$, yielding
$d_{i,i-1}=0$;

    3. (c)

$j=i+1$; Then $b_{i-2,j}=b_{j,i-2}=(2i^{3}-3i^{2}-11i+18)/6$,
$b_{i-1,j}=b_{j,i-1}=(2i^{3}-5i+3)/6$, $b_{i,j}=b_{j,i}=(2i^{3}+3i^{2}+i)/6$,
$b_{i+1,j}=(2i^{3}+6i^{2}+7i+3)/6$, $b_{i+2,j}=(2i^{3}+9i^{2}+13i+6)/6$,
yielding $d_{i,i+1}=0$.

    4. (d)

$j\geq i-1$; Then $b_{i-2,j}=b_{j,i-2}=(3j(i-2)^{2}+(i-2)-(i-2)^{3})/6$,
$b_{i-1,j}=b_{j,i-1}=(3j(i-1)^{2}+(i-1)-(i-1)^{3})/6$,
$b_{i,j}=b_{j,i}=\frac{3ji^{2}+i-i^{3}}{6}$,
$b_{i+1,j}=(3j(i+1)^{2}+(i+1)-(i+1)^{3})/6$,
$b_{i+2,j}=(3j(i+2)^{2}+(i+2)-(i+2)^{3})/6$, yielding $d_{i,j}=0$.

  2. (ii)

The case $i=1$.

For $j=1$, $b_{1,1}=3/6$, $b_{2,1}=1$, and $b_{3,1}=3/2$; Thus,
$d_{i,j}=d_{1,1}=7b_{1,1}-4b_{2,1}+b_{3,1}=1$.

For $j>1$, we have $b_{1,j}=b_{j,1}=j/2$, $b_{2,j}=b_{j,2}=2j-1$, and
$b_{3,j}=b_{j,3}=\frac{9j}{2}-4$; Thus, $d_{i,j}=7b_{1,j}-4b_{2,j}+b_{3,j}=0$.

  3. (iii)

The case $i=2$, with ${d_{i,j}}=-4b_{1,j}=6b_{2,j}-4b_{3,j}+b_{4,j}$.

For $j=2$, we have $b_{1,2}=b_{2,1}=1$, $b_{2,2}=3$, $b_{3,2}=5$, $b_{4,2}=7$;
Thus, $d_{2,2}=1$.

For $j\neq 2$, then $b_{1,j}=j/2$, $b_{2,j}=2j-1$, $b_{3,j}=\frac{9j}{2}-4$,
and $b_{4,j}=8j-10$. We have $d_{i,j}=0$.

  4. (iv)

For the case $i\in\\{n-1,n\\}$, similar computations using (14) complete the
proof.

∎

From now on, we shall use $a_{i,j}^{-1}$ to denote the $(i,j)$-entry of
$A^{-1}$, the inverse of $A$; thus, $a^{-1}_{i,j}=b_{i,j}$.

The following corollary is a consequence of Theorem 2.1.

######  Corollary 2.2.

The inverse of $A$ is a positive matrix; i.e., $A^{-1}>0$, implying
$a^{-1}_{i,j}>0$.

###### Proof.

By Theorem 2.1 it follows that $a^{-1}_{n,n}=n(2n^{2}+1)/12$ is positive.
Notice that, for $i\geq j$,
$a_{i,j}^{-1}=\displaystyle\frac{3ij^{2}+j-j^{3}}{6}\geq\frac{3j^{3}+j-j^{3}}{6}>0$.
Consequently, entries determined by the other 2 parts of Theorem 2.1 are also
positive. ∎

The above positivity result is important in the context of the fixed-point
iteration we devise to solve the nonlinear system (5). Consider the iteration

$\displaystyle\bm{u}^{\ell}=h^{4}KA^{-1}\exp(-\bm{u}^{\ell-1}),\ell=1,2.\dots$
(15)

Since $A^{-1}>0$ (Corollary 2.2), the recipe (15) generates a sequence of
positive vectors $\\{\bm{u}^{\ell}\\}$, if started with $\bm{u}^{0}>\bm{0}$.
As the solution of this type of boundary-value problem is a nonnegative
function (c.f., Section 1; see also later for the finite-difference equation
case), if the above iteration converges, it converges to a positive solution.

Let $p\in\\{1,2,\infty\\}$. Our starting point for the convergence analysis is
the relation, with $\bm{u}^{0}>\bm{0}$,

$\displaystyle\|\bm{u}^{\ell}-\bm{u}^{\ell-1}\|_{p}$
$\displaystyle=\|h^{4}KA^{-1}(\exp(-\bm{u}^{\ell-1})-\exp(-\bm{u}^{\ell-2}))\|_{p}$
$\displaystyle=h^{4}K\|A^{-1}\left(\exp(-\bm{u}^{\ell-2})+G(\bm{u}^{\ell-1}-\bm{u}^{\ell-2})-\exp(-\bm{u}^{\ell-2})\right)\|_{p}$
$\displaystyle=h^{4}K\|A^{-1}G(\bm{u}^{\ell-1}-\bm{u}^{\ell-2})\|_{p},$

where $G=-\text{diag}(\exp(-\xi_{1}),\dots,\exp(-\xi_{n}))$, such that the
vector
$\bm{\xi}=[\xi_{i}]_{i=1,n}\in\mathcal{B}=\\{\bm{x}\in\mathbb{R}^{n}:\|\bm{x}-\bm{u}^{\ell-2}\|_{p}<\|\bm{u}^{\ell-1}-\bm{u}^{\ell-2}\|_{p}\\}$.
Since $\\{\bm{u}^{\ell}\\}$ is a sequence of positive vectors, $\bm{\xi}$ is
also a positive vector, and consequently the diagonal entries of $G$ are
strictly less than 1. Thus, $\|G\|_{p}<1$, and

$\displaystyle\|\bm{u}^{\ell}-\bm{u}^{\ell-1}\|_{p}$ $\displaystyle\leq$
$\displaystyle
h^{4}K\|A^{-1}\|_{p}\|G\|_{p}(\bm{u}^{\ell-1}-\bm{u}^{\ell-2})\|_{p}$ (16)
$\displaystyle<$ $\displaystyle
h^{4}K\|A^{-1}\|_{p}\|\bm{u}^{\ell-1}-\bm{u}^{\ell-2}\|_{p}.$

We define $L_{p}$ to be

$L_{p}=h^{4}K\|A^{-1}\|_{p}.$ (17)

Convergence guarantee of the fixed point iteration (15) requires $L_{p}<1$,
which in turn, for given $K$ and chosen $h$, requires that

$\|A^{-1}\|_{p}<1/(h^{4}K).$ (18)

######  Lemma 2.3.

For the inverse of $A$ in Theorem 2.1, the following holds true:

$a_{i_{1},j}^{-1}>a_{i_{2},j}^{-1},\quad\forall i_{1}>i_{2}>j,\text{ with
}i_{1},i_{2},j\in\\{1,2,\dots,n\\}.$

###### Proof.

From Theorem 2.1 it follows that $a^{-1}_{i_{1},j}=(3i_{1}j^{2}+j-j^{3})/6$
and $a_{i_{2},j}^{-1}=(3i_{2}j^{2}+j-j^{3})/6$, thus, one can notice that
$a_{i_{1},j}^{-1}>a_{i_{2},j}^{-1}$, for $i_{1}>i_{2}>j$. ∎

######  Theorem 2.4.

Let $A\in\mathbb{R}^{n\times n}$ be given in (13), with $n\geq 5$. Then

$\|A^{-1}\|_{p}=\begin{cases*}(n^{4}-n^{2})/8,\quad\text{if }p=1,\\\\[5.0pt]
(n^{4}+n^{2})/8,\quad\text{if }p=\infty.\end{cases*}$

###### Proof.

For $p=1$ case, it follows from Lemma 2.3 that

$\|A^{-1}\|_{1}=\max_{1\leq j\leq
n}\sum_{i=1}^{n}|a_{i,j}^{-1}|=\max\left\\{\sum_{i=1}^{n}|a^{-1}_{i,n-1}|,\sum_{i=1}^{n}|a^{-1}_{i,n}|\right\\}.$

We have

$\sum_{i=1}^{n}|a_{i,n-1}^{-1}|=\sum_{i=1}^{n}|a_{n-1,i}^{-1}|=\sum_{i=1}^{n}\frac{3(n-1)i^{2}+i-i^{3}}{6}=\frac{n^{4}-n^{2}}{8}.$

We can now proceed similarly:

$\sum_{i=1}^{n}|a_{i,n}^{-1}|=\frac{1}{2}\sum_{i=1}^{n}|a_{n,i}^{-1}|=\frac{1}{2}\sum_{i=1}^{n}\frac{3ni^{2}+i-i^{3}}{6}=\frac{3n^{4}+4n^{3}+3n^{2}+2n}{48}.$

From the above results,

$\sum_{i=1}^{n}|a_{i,n-1}^{-1}|-\sum_{i=1}^{n}|a_{i,n}^{-1}|=\frac{3n^{4}-4n^{3}-9n^{2}-2n}{48}>0$

for $n\geq 5$. Therefore,

$\max_{1\leq j\leq
n}\sum_{i=1}^{n}|a_{i,j}^{-1}|\leq\sum_{i=1}^{n}|a_{i,n-1}^{-1}|=\frac{n^{4}-n^{2}}{8}=\|A^{-1}\|_{1}.$

Next for $p=\infty$ using the Lemma 2.3,

$\displaystyle\|A^{-1}\|_{\infty}$ $\displaystyle=$ $\displaystyle\max_{1\leq
i\leq
n}\sum_{j=1}^{n}|a_{i,j}^{-1}|=\sum_{j=1}^{n}|a_{n,j}^{-1}|=\sum_{j=1}^{n-1}|a_{n,j}^{-1}|+a_{n,n}^{-1}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n-1}\frac{3nj^{2}+j-j^{3}}{6}+\frac{n(2n^{2}+1)}{12}$
$\displaystyle=$ $\displaystyle\frac{n^{4}+n^{2}}{8}.$

∎

Using Hölder’s inequality,

$\|A^{-1}\|_{2}\leq\sqrt{\|A^{-1}\|_{1}\|A^{-1}\|_{\infty}}=\frac{1}{8}\sqrt{n^{8}-n^{4}}\leq\frac{1}{8}n^{4}.$

We conclude this section by the characterization of the finite difference
solution of the system (5). Because $A\bm{u}=h^{4}K\exp(-\bm{u})>\bm{0}$, for
the last row of the system,

$\displaystyle 2u_{n-2}-4u_{n-1}+2u_{n}>0\quad\Longrightarrow\quad
u_{n}-u_{n-1}>u_{n-1}-u_{n-2}.$ (19)

From the $(n-1)$-th row, with $u_{n-3}-4u_{n-2}+5u_{n-1}-2u_{n}>0$, we have

$\displaystyle
u_{n-1}-u_{n-2}>u_{n-2}-u_{n-3}+2u_{n-2}-4u_{n-1}+2u_{n}>u_{n-2}-u_{n-3},$
(20)

after using the inequality (19). Furthermore, this row leads to

$4(u_{n-1}-u_{n-2})>u_{n}-u_{n-3}+u_{n}-u_{n-1}>u_{n}-u_{n-3}+u_{n-1}-u_{n-2},$

after again using (19), which in turn yields

$\displaystyle 3(u_{n-1}-u_{n-2})>u_{n}-u_{n-3}.$ (21)

We then have the following lemma:

######  Lemma 2.5.

For the inequality $A\bm{u}>\mathbf{0}$, with $A$ given by (13), the following
inequalities hold, with $j=i+2$ and $i=3\dots,n-1$ the rows of $A$:

$u_{j+2}-u_{j+1}>u_{j+1}-u_{j},\quad 3(u_{j+2}-u_{j+1})>u_{j+3}-u_{j}.$

###### Proof.

We have proved the inequalities for $j=n-3$, which comes from the $(n-1)$-th
row of $A\bm{u}>\mathbf{0}$. Now suppose that they hold also for
$j=n-3,n-4,\dots,k+1$. Associated with $j=k$ is the inequality
$u_{k}-4u_{k+1}+6u_{k+2}-4u_{k+3}+u_{k+4}>0$ from the $(k+2)$-th row of
$A\bm{u}>\mathbf{0}$, which gives

$\displaystyle u_{k+2}-u_{k+1}$ $\displaystyle>$ $\displaystyle
3u_{k+1}-u_{k}-5u_{k+2}+4u_{k+3}-u_{k+4}$ $\displaystyle=$ $\displaystyle
u_{k+1}-u_{k}+\left[4(u_{k+3}-u_{k+2})+u_{k+1}-u_{k+4}+u_{k+1}-u_{k+2}\right]$
$\displaystyle>$ $\displaystyle
u_{k+1}-u_{k}+\left[3(u_{k+3}-u_{k+2})+u_{k+1}-u_{k+4}\right]$
$\displaystyle>$ $\displaystyle u_{k+1}-u_{k}$

by assumption. Next, note that
$u_{k}-4u_{k+1}+6u_{k+2}-4u_{k+3}+u_{k+4}=3(u_{k+2}-u_{k+1})+u_{k}-u_{k+3}-[3(u_{k+3}-u_{k+2})+u_{k+1}-u_{k+4}]>0$.
Thus, $3(u_{k+2}-u_{k+1})+u_{k}-u_{k+3}>3(u_{k+3}-u_{k+2})+u_{k+1}-u_{k+4}>0$,
by assumption. ∎

######  Theorem 2.6.

The solution of the finite difference system (5) is a nonnegative vector
$\bm{u}$, with increasing $u_{i}$.

###### Proof.

On the nodes $i=0,1$, approximation to the differential term leads to

$u_{i-2}-4u_{i-1}+6u_{i}-4u_{i+1}+u_{i+2}>0,$

which is of the same structure as the $i=3,\dots,n-2$ rows of $A$. By Lemma
2.5,

$u_{2}-u_{1}>u_{1}-u_{0},\quad u_{1}-u_{0}>u_{0}-u_{-1}.$

Therefore,

$u_{n}-u_{n-1}>u_{n-1}-u_{n-2}>\cdots>u_{2}-u_{1}>u_{1}-u_{0}>u_{0}-u_{-1}.$

With $u_{0}=0$ (from the boundary condition $\phi(0)=\phi_{0}=0)$ and
$u_{-1}=u_{1}$ (from using central finite differencing on
$\phi^{\prime}(0)=0$), from the most right inequality, we get $u_{1}>0=u_{0}$.
Also, $u_{2}-u_{1}>u_{1}-u_{0}>0$; thus $u_{2}>u_{1}$. In general, we have
$u_{i+1}>u_{i}$, $i=1,\dots,n-1$. ∎

##  3 The case with clamped-clamped boundary conditions

In this section, we consider the case with the boundary conditions (3).
Conditions at $x=1$ are treated in the same way as at $x=0$, leading to (5),
but now with $\bm{u}=(u_{1},\dots,u_{n-1})^{T}\in\mathbb{R}^{n-1}$ and
$A\in\mathbb{R}^{(n-1)\times(n-1)}$ given by

$\displaystyle A=\left[\begin{array}[]{rrrrrrr}7&-4&1&0&&\cdots&0\\\
-4&6&-4&\ddots&\ddots&&\vdots\\\ 1&-4&\ddots&\ddots&&\ddots&\\\
0&\ddots&\ddots&&\ddots&\ddots&0\\\ &\ddots&&&\ddots&-4&1\\\
\vdots&&\ddots&\ddots&-4&6&-4\\\ 0&\cdots&&0&1&-4&7\end{array}\right].$ (29)

However, to simplify our notation, we shall consider the case where
$\bm{u}\in\mathbb{R}^{n}$ and $A\in\mathbb{R}^{n\times n}$ in the subsequent
analysis; in this case, $h=1/(n+1)$.

In constrast to (13), the matrix (29) is centrosymmetric and near Toeplitz.
Furthermore, it admits the rank-2 decomposition as follows:

$\displaystyle A=T^{2}+UU^{t},$ (30)

where $T=\text{tridiag}_{n}(-1,2,-1)$ is an $n\times n$ tridiagonal symmetric
Toeplitz matrix, and

$\displaystyle U=\left[\begin{array}[]{cc}\sqrt{2}&0\\\ 0&0\\\
\vdots&\vdots\\\ 0&\sqrt{2}\end{array}\right]\in\mathbb{R}^{n\times 2}.$ (35)

$T$ is a symmetric M-matrix, with positive inverse given explicitly by (see,
e.g., [10])

$\displaystyle[T^{-1}]_{ij}=\begin{cases}\frac{j}{n+1}(n-(i-1)),&i\geq
j,\\\\[5.0pt] \frac{i}{n+1}(n-(j-1)),&i<j.\end{cases}$ (36)

$A$ is symmetric positive definite because $T^{2}=T^{T}T$ (and $UU^{t}$) is
symmetric positive (semi) definite. The inverse of $A$ can be computed by
applying the Sherman-Morrison formula on (30):

$\displaystyle A^{-1}$
$\displaystyle=T^{-2}-T^{-2}U(I_{2}+U^{t}T^{-2}U)^{-1}U^{t}T^{-2}$
$\displaystyle=T^{-1}(I-T^{-1}U(I_{2}+U^{t}T^{-2}U)^{-1}U^{t}T^{-1})T^{-1}$
$\displaystyle=T^{-t}(I-T^{-1}U(I_{2}+U^{t}T^{-2}U)^{-1}(T^{-1}U)^{t})T^{-1}.$
(37)

Because $A^{-1}$ is symmetric positive definite, the middle term on the right-
hand side $I-T^{-1}U(I_{2}+U^{T}T^{-2}U)(T^{-1}U)^{T}$ is also symmetric
positive definite. Rewriting (30) as

$A=T^{2}+UU^{t}=T^{t}(I+T^{-1}U(T^{-1}U)^{t})T,$

clearly

$(I+T^{-1}U(T^{-1}U)^{t})^{-1}=I-T^{-1}U(I_{2}+U^{t}T^{-2}U)^{-1}(T^{-1}U)^{t}=:M.$

Note that, with (35) and (36),

$\displaystyle T^{-1}U=\frac{\sqrt{2}}{n+1}\left[\begin{array}[]{cc}n&1\\\
n-1&2\\\ \vdots&\vdots\\\ 2&n-1\\\ 1&n\end{array}\right].$ (43)

Direct computation yields

$\displaystyle I_{2}+U^{t}T^{-2}U$ $\displaystyle=$
$\displaystyle\frac{2}{(n+1)^{2}}\left[\begin{array}[]{ccc}\displaystyle\frac{(n+1)^{2}}{2}+\sum_{k=1}^{n}k^{2}&&\displaystyle\sum_{k=1}^{n}(n-(k-1))k\\\
&&\\\
\displaystyle\sum_{k=1}^{n}(n-(k-1))k&&\displaystyle\frac{(n+1)^{2}}{2}+\sum_{k=1}^{n}k^{2}\end{array}\right]$
(47) $\displaystyle=$
$\displaystyle\frac{1}{\gamma}\left[\begin{array}[]{ccc}\displaystyle\gamma+\tau&&\displaystyle\gamma
n-\tau\\\ &&\\\ \displaystyle\gamma
n-\tau&&\displaystyle\gamma+\tau\end{array}\right],$ (51)

where $\tau=\frac{2n^{3}+3n^{2}+n}{6}$ and $\gamma=\frac{(n+1)^{2}}{2}$. Its
inverse is given by

$\displaystyle(I_{2}+U^{T}T^{-2}U)^{-1}=\frac{1}{\delta}\left[\begin{array}[]{ccc}\displaystyle\gamma+\tau&&\displaystyle-\gamma
n+\tau\\\ &&\\\ \displaystyle-\gamma
n+\tau&&\displaystyle\gamma+\tau\end{array}\right],$ (55)

where $\det(I_{2}+U^{T}T^{-2}U)=\frac{1}{3}(n^{2}+2n+3)>0$ and
$\delta=(n+1)(2\tau+\gamma(1-n))=\frac{1}{6}(n+1)^{2}(n^{2}+2n+3)$.

Let $M=[m_{ij}]_{i,j=1,n}$. Using (43) and (55), we have, for $i\neq j$,

$\displaystyle m_{ij}$ $\displaystyle=$
$\displaystyle-\frac{2}{\delta(n+1)^{2}}[n-(i-1)\quad
i]\left[\begin{array}[]{ccc}\displaystyle\gamma+\tau&&\displaystyle-\gamma
n+\tau\\\ &&\\\ \displaystyle-\gamma
n+\tau&&\displaystyle\gamma+\tau\end{array}\right]\left[\begin{array}[]{c}n-(j-1)\\\
j\end{array}\right]$ $\displaystyle=$ $\displaystyle
q_{0}(n)+q_{1}(n)(i+j)+q_{2}(n)ij,$

where

$\displaystyle q_{0}(n)$ $\displaystyle=$
$\displaystyle-\frac{4n^{2}+8n+6}{(n+1)(n^{2}+2n+3)},$ $\displaystyle
q_{1}(n)$ $\displaystyle=$ $\displaystyle\frac{6}{n^{2}+2n+3},$ $\displaystyle
q_{2}(n)$ $\displaystyle=$ $\displaystyle-\frac{12}{(n+1)(n^{2}+2n+3)}.$

One can verify that $m_{ij}$ change signs. Thus $M$ is not an M-matrix.

For $i=j$,

$\displaystyle m_{ii}$ $\displaystyle=$ $\displaystyle
1-\frac{2}{\delta(n+1)^{2}}[n-(i-1)\quad
i]\left[\begin{array}[]{ccc}\displaystyle\gamma+\tau&&\displaystyle-\gamma
n+\tau\\\ &&\\\ \displaystyle-\gamma
n+\tau&&\displaystyle\gamma+\tau\end{array}\right]\left[\begin{array}[]{c}n-(i-1)\\\
i\end{array}\right]$ $\displaystyle=$
$\displaystyle\frac{n^{3}-n^{2}-3n-3}{(n+1)(n^{2}+2n+3)}+\frac{12}{n^{2}+2n+3}i-\frac{12}{(n+1)(n^{2}+2n+3)}i^{2}$
$\displaystyle>$ $\displaystyle 0,$

for $n\geq 1$.

######  Theorem 3.1.

The inverse of $A$ given by (29) is a positive matrix. Furthermore, let
$\alpha=n+1-i$, $\beta=j\alpha/(6(n+1)(n^{2}+2n+3))$, and
$\varepsilon=3(1+\alpha(n+1))(1+(i-j)j)$. The entries of $A^{-1}$ are

  * •

$a^{-1}_{ij}=\beta(\varepsilon+(j^{2}-1)(2\alpha^{2}+1))$, for $i\geq j$

  * •

$a^{-1}_{ij}=a^{-1}_{ji}$, otherwise.

###### Proof.

Let $A^{-1}=[a^{-1}_{ij}]$, with $A^{-1}=T^{-1}MT^{-1}$. Denote by
$\mathbf{y}_{j}=[y_{k,j}]_{k=1,n}=MT_{j}^{-1}$, the product of $M$ and the
$j$-th column of $T^{-1}$. For $i\geq j$,

$\displaystyle\mathbf{y}_{j}=\frac{1}{n+1}\begin{bmatrix}(n+1-j)(m_{1,1}+2m_{1,2}+\cdots+jm_{1,j})+j(m_{1,j+1}(n-j)+\cdots+m_{1,n})\\\
\vdots\\\
(n+1-j)(m_{j,1}+2m_{j,2}+\cdots+jm_{j,j})+j(m_{j,j+1}(n-j)+\cdots+m_{j,n})\\\
\vdots\\\
(n+1-j)(m_{n,1}+2m_{n,2}+\cdots+jm_{n,j})+j(m_{n,j+1}(n-j)+\cdots+m_{n,n})\end{bmatrix}.$

Using $m_{i,j}$, for $i\leq j$, we have, with
$m^{*}_{i,j}=q_{0}+q_{1}(i+j)+q_{2}ij$,

$\displaystyle y_{i,j}$ $\displaystyle=$ $\displaystyle\frac{(n-(j-1))i}{n+1}$
$\displaystyle+$
$\displaystyle\frac{n+1-j}{n+1}(m^{*}_{i,1}+2m^{*}_{i,2}+\cdots+jm^{*}_{i,j})+\frac{j}{n+1}(m^{*}_{i,j+1}(n-j)+\cdots+m^{*}_{i,n})$
$\displaystyle=$ $\displaystyle\frac{(n-(j-1))i}{n+1}$ $\displaystyle+$
$\displaystyle\frac{q_{0}+q_{1}i}{n+1}((n+1-j)(1+\cdots+j)+j(n-j+\cdots+1))$
$\displaystyle+$
$\displaystyle\frac{q_{1}+q_{2}i}{n+1}((n+1-j)(1^{2}+\cdots+j^{2})+j((n-j)(j+1)+(n-j-1)(j+2)+\cdots+n))$
$\displaystyle=$
$\displaystyle\frac{1}{n+1}\left((n-(j-1))i+r_{0}+r_{1}i\right),$

where

$r_{0}=-\frac{j(n+1-j))((n+1)(n+1-j)+1)}{n^{2}+2n+3}$

and

$r_{1}=\frac{j(n+1-j)(n+1-2j)}{n^{2}+2n+3}.$

Using similar calculation for $i>j$, we get

$y_{i,j}=\frac{1}{n+1}\begin{cases}\displaystyle
r_{0}+r_{1}i+(n+1-j)i=r_{0}+\left(r_{1}-j\right)i+(n+1)i,&i\leq j,\\\
\displaystyle
r_{0}+r_{1}i+(n+1-i)j=r_{0}+\left(r_{1}-i\right)i+(n+1)j,&i>j;\end{cases}$

hence,

$\mathbf{y}_{j}=\frac{1}{n+1}\left(r_{0}\begin{bmatrix}1\\\ \vdots\\\ 1\\\
1\\\ \vdots\\\ 1\end{bmatrix}+\left(r_{1}-j\right)\begin{bmatrix}1\\\
\vdots\\\ j-1\\\ j\\\ \vdots\\\ n\end{bmatrix}+(n+1)\begin{bmatrix}1\\\
\vdots\\\ j-1\\\ j\\\ \vdots\\\ j\end{bmatrix}\right).$

Consider the $i$-th row of $T^{-1}$:

$T_{i}^{-t}=\frac{1}{n+1}\begin{bmatrix}n+1-i&2(n+1-i)&\cdots&i(n+1-i)&i(n-i)&\cdots&i\end{bmatrix}.$

We have

$\displaystyle a^{-1}_{ij}$ $\displaystyle=$ $\displaystyle
T_{i}^{-t}MT_{j}^{-1}=T_{i}^{-t}\mathbf{y}_{j}$ $\displaystyle=$
$\displaystyle
r_{0}\frac{i(n+1-i)}{2}+\left(r_{1}-\frac{j}{n+1}\right)\frac{i(n+1-i)(n+1+i)}{6}$
$\displaystyle+$ $\displaystyle\frac{j(n+1-i)(3i(n+1)+1-j^{2})}{6(n+1)}$
$\displaystyle=$ $\displaystyle\beta(\varepsilon+(j^{2}-1)(2\alpha^{2}+1)),$

where

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle\displaystyle{n+1-i},$
$\displaystyle\beta$ $\displaystyle=$
$\displaystyle\displaystyle\frac{j(n+1-i)}{6(n+1)(n^{2}+2n+3)},$
$\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle\displaystyle{3(1+\alpha+n\alpha)(1+ij-j^{2})}.$

Notice that $\alpha,\beta,\varepsilon>0$, $\forall i,j=1,\dots,n$. With $i\geq
j$ and $j^{2}>j^{2}-1$,

$\displaystyle a_{ij}^{-1}$ $\displaystyle>$
$\displaystyle\beta(j^{2}-1)(2\alpha^{2}+1)$ $\displaystyle\geq$
$\displaystyle 0.$

∎

By Theorem 3.1, starting from $\bm{u}^{0}>\bm{0}$, the fixed-point iteration
(15) is guaranteed to generate a sequence of positive vectors.

In the sequel, we present two ways of constructing an estimate for norms of
the inverse of $A$. The first approach is based on the factorization
$A^{-1}=T^{-t}M^{-1}T^{-1}$ in (37). The result is presented in the next
theorem.

######  Theorem 3.2.

For $p\in\\{1,2,\infty\\}$,

$\|A^{-1}\|_{p}\leq(n+1)^{4}/32.$

###### Proof.

$\|A^{-1}\|_{p}\leq\|T^{-1}\|_{p}\|M\|_{p}\|T^{-1}\|_{p}=\|T^{-1}\|_{p}^{2}\|M\|_{p}.$

Note that $\|T^{-1}\|_{1}=\|T^{-1}\|_{\infty}$, due to symmetry. Thus, we
shall consider only $\|T^{-1}\|_{1}$. Using (36),

$\displaystyle\sum_{j=1}^{n}|T^{-1}_{ij}|$ $\displaystyle=$
$\displaystyle\frac{1}{n+1}\left[(n-(i-1))\sum_{j=1}^{i-1}j+i\sum_{j=1}^{n-(i-1)}j\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2(n+1)}\left[(n+1)^{2}i-(n+1)i^{2}\right].$

The maximum of the rowsum is then attained for $i=(n+1)/2$. Thus,

$\displaystyle\|T^{-1}\|_{1}=\max_{1\leq i\leq
n}\sum_{j=1}^{n}|T^{-1}_{ij}|\leq\frac{(n+1)^{2}}{8},$ (68)

with equality holding when $n$ is odd.

We now estimate the 1-norm of $M$. Let
$\widetilde{m}_{ij}=q_{0}(n)+q_{1}(n)(i+j)+q_{2}(n)ij$, $\forall
i,j=1,\dots,n$ and consider $\sum_{j=1}^{n}|\widetilde{m}_{ij}|$. For a fixed
$i$, $\widetilde{m}_{ij}$ can be viewed as a linear function of $j$.
$\sum_{j=1}^{n}|\widetilde{m}_{ij}|$ can then be viewed as the rectangular
rules that approximate the area made by the function $\widetilde{m}_{ij}$ and
the $j$-axis. In this case, treating $j\in[0,n+1]\subset\mathbb{R}$,

$\sum_{j=1}^{n}|\widetilde{m}_{ij}|\leq\int_{j=0}^{n+1}|\widetilde{m}_{ij}|dj=\frac{1}{2}(|\widetilde{m}_{i,0}|+|\widetilde{m}_{i,n+1}|)(n+1),$

where $\widetilde{m}_{i,0}=-(4n^{2}+6n(1-i)+8-6i)/[(n+1)(n^{2}+2n+3)]$ and
$\widetilde{m}_{i,n+1}=(2n^{2}+6n-2-6i(n+1))/[(n+1)(n^{2}+2n+3)]$.

Since the matrix $\widetilde{M}=[\widetilde{m}_{ij}]$ is persymmetric, we just
need to consider $i=1,\dots,(n+1)/2$. Then,

$\displaystyle\sum_{j=1}^{n}|\widetilde{m}_{ij}|$ $\displaystyle\leq$
$\displaystyle\max_{i}\frac{1}{2}(|\widetilde{m}_{i,0}|+|\widetilde{m}_{i,n+1}|)(n+1)=\frac{1}{2}\frac{2(n^{2}+5)}{(n+1)(n^{2}+2n+3)}(n+1)$
$\displaystyle=$ $\displaystyle\frac{n^{2}+5}{n^{2}+2n+3}.$

Now,

$\displaystyle\sum_{j=1}^{n}|m_{ij}|$ $\displaystyle=$
$\displaystyle\sum_{j=1,j\neq 1}^{n}|m_{ij}|+|m_{ii}|=\sum_{j=1,j\neq
1}^{n}|\widetilde{m}_{ij}|+|1+\widetilde{m}_{ii}|$ $\displaystyle\leq$
$\displaystyle 1+|\widetilde{m}_{ii}|+\sum_{j=1,j\neq
i}^{n}|\widetilde{m}_{ij}|$ $\displaystyle=$ $\displaystyle
1+\sum_{j=1}^{n}|\widetilde{m}_{ij}|.$

Thus, for $n\geq 1$,

$\displaystyle\|M\|_{1}=\max_{i}\sum_{j=1}^{n}|m_{ij}|$ $\displaystyle\leq$
$\displaystyle 1+\max_{i}\sum_{j=1}^{n}|\widetilde{m}_{ij}|$
$\displaystyle\leq$ $\displaystyle 1+\frac{n^{2}+5}{n^{2}+2n+3}$
$\displaystyle\leq$ $\displaystyle 2,$

since $n^{2}+5<n^{2}+2n+3$ for $n\geq 1$.

Combining with $\|T^{-1}\|_{1}$, we get the desired result. Furthermore, using
Hölder’s inequality,
$\|A^{-1}\|_{2}\leq\sqrt{\|A^{-1}\|_{1}\|A^{-1}\|_{\infty}}\leq(n+1)^{4}/32$.
∎

The second approach uses the knowledge of the entries of $A^{-1}$ in Theorem
3.1. Tedious calculation results in exact norms in some cases, and hence much
stronger estimates than the previous estimates.

######  Theorem 3.3.

For $p\in\\{1,2,\infty\\}$,

$\|A^{-1}\|_{p}\leq(n+1)^{2}\left((n+1)^{2}+8\right)/384.$

If $n$ is odd, then the equality holds for $p\in\\{1,\infty\\}$.

###### Proof.

We shall first consider the case $p=\infty$. In this case, by using
$a_{i,j}^{-1}>0$,

$\displaystyle\|A^{-1}\|_{\infty}$ $\displaystyle=$
$\displaystyle\max_{i}\sum_{j=1}^{n}|a^{-1}_{i,j}|=\max_{i}\sum_{j=1}^{n}a^{-1}_{i,j}$

For $i=1,\dots,n$,

$\displaystyle\sum_{j=1}^{n}a_{i,j}^{-1}=\sum_{j=1}^{i}a^{-1}_{i,j}+\sum_{j=i+1}^{n}a^{-1}_{i,j}=\sum_{j=1}^{i}a^{-1}_{i,j}+\sum_{k=1}^{n-i}a^{-1}_{k,i},$

because of the centrosymmetry of $A^{-1}$. Calculating each sum using the
formula for the entries $a^{-1}_{ij}$, we get

$\displaystyle\sum_{j=1}^{i}a^{-1}_{i,j}$ $\displaystyle=$
$\displaystyle\widehat{\delta}^{-1}\left[C^{i}_{1}\sum_{j=1}^{i}j+C^{i}_{2}\sum_{j=1}^{i}j^{2}+C^{i}_{3}\sum_{j=1}^{i}j^{3}\right]$
$\displaystyle=$
$\displaystyle\widehat{\delta}^{-1}\left[C^{i}_{1}\frac{i^{2}+i}{2}+C^{i}_{2}\frac{2i^{3}+3i^{2}+i}{6}+C^{i}_{3}\frac{i^{4}+2i^{3}+i^{2}}{4}\right],$

where $\widehat{\delta}=6(n+1)(n^{2}+2n+3)$ and

$\displaystyle C^{i}_{1}$ $\displaystyle=$ $\displaystyle
n^{3}+3n^{2}-3i^{2}n+5n+2i^{3}-3i^{2}-2i+3,$ $\displaystyle C^{i}_{2}$
$\displaystyle=$ $\displaystyle
3in^{3}-6i^{2}n^{2}+9in^{2}+3i^{3}n-12i^{2}n+12in+3i^{3}-9i^{2}+6i,$
$\displaystyle C^{i}_{3}$ $\displaystyle=$
$\displaystyle-n^{3}-3n^{2}+3i^{2}n-5n-2i^{3}+3i^{2}+2i-3.$

Also,

$\displaystyle\sum_{k=1}^{n-i}a^{-1}_{k,i}$ $\displaystyle=$
$\displaystyle\widehat{\delta}^{-1}\left[C^{k}_{1}\sum_{k=1}^{n-i}k+C^{k}_{2}\sum_{k=1}^{n-i}k^{2}+C^{k}_{3}\sum_{k=1}^{n-i}k^{3}\right]$
$\displaystyle=$
$\displaystyle\widehat{\delta}^{-1}\left[C^{k}_{1}\frac{(n-i)^{2}+(n-i)}{2}+C^{k}_{2}\frac{2(n-i)^{3}+3(n-i)^{2}+(n-i)}{6}\right]$
$\displaystyle+$
$\displaystyle\widehat{\delta}^{-1}C^{k}_{3}\frac{(n-i)^{4}+2(n-i)^{3}+(n-i)^{2}}{4},$

where

$\displaystyle C^{k}_{1}$ $\displaystyle=$ $\displaystyle
3i^{2}n-2i^{3}+3i^{2}+2i,$ $\displaystyle C^{k}_{2}$ $\displaystyle=$
$\displaystyle 3i^{2}n^{2}-3i^{3}n+6i^{2}n+3in-3i^{3}+3i,$ $\displaystyle
C^{k}_{3}$ $\displaystyle=$ $\displaystyle-3i^{2}n+2i^{3}-3i^{2}-2i.$

Assuming that $i\in[1,n]\subset\mathbb{R}$, the maximum of the rowsum is
obtained from the condition
$\displaystyle\frac{d}{di}\sum_{j=1}^{n}a_{i,j}^{-1}=0$. In this regard, we
have

$\displaystyle\frac{d}{di}\sum_{j=1}^{n}a_{i,j}^{-1}=\widehat{\delta}^{-1}\left[C^{\prime}_{0}+C^{\prime}_{1}i+C^{\prime}_{2}i^{2}+C^{\prime}_{3}i^{3}\right]=0,$

where

$\displaystyle C^{\prime}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}n^{4}+2n^{3}+4n^{2}+4n+\frac{3}{2},$ $\displaystyle
C^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}n^{5}+\frac{5}{2}n^{4}+5n^{3}+5n^{2}+\frac{1}{2}n-\frac{3}{2},$
$\displaystyle C^{\prime}_{2}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}n^{4}-6n^{3}-12n^{2}-12n-\frac{9}{2},$
$\displaystyle C^{\prime}_{3}$ $\displaystyle=$ $\displaystyle
n^{3}+3n^{2}+5n+3.$

The only acceptable solution of the above equation is $i=(n+1)/2$. The other
solutions are rejected: $i=-\frac{1}{2}(\sqrt{n^{2}+2n+5}-(n+1))<0$ and
$i=\frac{1}{2}(\sqrt{n^{2}+2n+5}+(n+1))>n+1>n$. One can verify that
$i=(n+1)/2$ maximizes the rowsum.

Let $n$ be odd. With $i=(n+1)/2$,

$\displaystyle\|A^{-1}\|_{\infty}$
$\displaystyle=\max_{i}\sum_{j=1}^{n}|a^{-1}_{i,j}|=\sum_{j=1}^{n}a^{-1}_{(n+1)/2,j}=\sum_{j=1}^{\frac{n+1}{2}}a^{-1}_{(n+1)/2,j}+\sum_{k=1}^{\frac{n-1}{2}}a^{-1}_{k,(n+1)/2}$
$\displaystyle=\left(n^{4}+4n^{3}+14n^{2}+20n+9\right)/384$
$\displaystyle=(n+1)^{2}((n+1)^{2}+8)/384.$

If $n$ is even, then $i=(n+1)/2$ is not a row of the matrix $A$; the maximum
of the rowsum will then be attained at $i=\lceil(n+1)/2\rceil$ or
$i=\lfloor(n+1)/2\rfloor$. Either case satisfies

$\displaystyle\|A^{-1}\|_{\infty}\leq(n+1)^{2}((n+1)^{2}+8)/384.$

Symmetry of $A^{-1}$ leads to $\|A^{-1}\|_{1}=\|A^{-1}\|_{\infty}$. Using
Hölder’s inequality, the above inequality holds also for $p=2$. ∎

Table 1 shows the computed norms of the inverse and compares them with the
estimate given by Theorem 3.3. For odd $n$ and $p\in\\{1,\infty\\}$ the norms
are exact. For even $n$, Theorem 3.3 gives an estimate that leads to a small
gap. This gap relative to the estimate becomes negligible with an increase in
$n$. To support this statement, the reader is referred to Fig. 1 and Fig. 2 in
log scales. The numerical tests are performed for all even $n$ from $10$ to
$1000$. The relative error is computed as
$|\|A^{-1}\|_{p}-UBound|/\|A^{-1}\|_{p}$, where
$UBound=(n+1)^{2}\left((n+1)^{2}+8\right)/384$ from Theorem 3.3. As shown in
Fig. 2 (left), the relative error decreases as $n$ increases for $p=1$ or
$p=\infty$. On the other hand, according to the numerical observation the
difference between $\|A\|_{2}$ and the upper bound become constant relative to
the norm as $n$ increases, see Fig. 2 (right).

Table 1: Computed $\|A^{-1}\|_{p}$ and the estimates, for the clamped-clamped case. $n$ | $p=$ | Upper bound from
---|---|---
| 1 | 2 | $\infty$ | Theorem 3.3
49 | 16,328 | 12,527 | 16,328 | 16,328
50 | 17,658 | 13,558 | 17,658 | 17,672
99 | 260,625 | 199,939 | 260,625 | 260,625
100 | 271,150 | 208,055 | 271,150 | 271,203
150 | 1,354,225 | 1,038,976 | 1,354,225 | 1,354,343

Figure 1: The upper bound and actual norm $p=1$ or $p=\infty$(left) and
$p=2$(right) in log scale

(a) $p=1$ or $p=\infty$(left)

(b) $p=2$(right)

Figure 2: The relative errors in log scale

For $n\geq 5$, the factor $(1+8/(n+1)^{2})/384\leq 11/3474$. So,
alternatively, if $K$ satisfies the condition in the above theorem, we can
have a simpler bound: $L_{p}<11/3474$. This factor approaches $1/384$ from
above as $n\to\infty$. Since the latter is slightly less than the former, for
a fixed $K$, one can expect a slight improvement of convergence by increasing
$n$.

##  4 Numerical Results

In this section, we present numerical results from solving (1) with (2) or (3)
using the fixed point method(15). We compare the observed convergence with the
theoretical bound given by (17) and Theorem 2.4 (for the clamped-free case) or
Theorem 3.3 (for the clamped-clamped case).The fixed point method (15) is
declared to have reached a convergence if
$\|\bm{u}^{\ell+1}-\bm{u}^{\ell}\|_{p}<10^{-6}$, where $p\in\\{1,2,\infty\\}$.
Solutions at convergence are shown in Figure 3 for the clamped-free and
clamped-clamped case, with $K=1$.

Figure 3: Solution at convergence with $K=1$, $n=100$: clamped-free (left),
clamped-clamped (right)

For both cases, the actual convergence rates are lower than the estimate
(Tables 2–5), with increasing gaps between the two as $K$ increases. As
$\|A^{-1}\|_{p}$ is exact, except for $p=2$, (due to the explicit inverse of
$A$), this suggests that the gap in the convergence rate is mainly due to the
estimate $\|G\|_{p}<1$. The numerical experiments suggest that the simple
fixed-point method (15) can be used for a wider range of $K$ than suggested by
the theoretical results. For instance, with $K=386$ and $n=99$, we have
$L_{p}=1.006$. The method still however converges to the solution at the
maximum rate of $0.5278$.

Table 2: Observed maximum convergence rate for clamped-free case, with $n=50$. In brackets are the theoretical rate based on Theorem 2.4. | $p=$
---|---
$K$ | $1$ |  | $2$ |  | $\infty$
$1/8$ | 0.010 | [0.016] |  | 0.010 | [0.016] |  | 0.010 | [0.017]
1 | 0.074 | [0.125] |  | 0.074 | [0.125] |  | 0.074 | [0.125]
8 | 0.400 | [1.000] |  | 0.400 | [1.000] |  | 0.402 | [1.000]
Table 3: Observed maximum convergence rate for clamped-free case, with $n=99$. In brackets are the theoretical rate based on Theorem 2.4. | $p=$
---|---
$K$ | $1$ |  | $2$ |  | $\infty$
$1/8$ | 0.010 | [0.016] |  | 0.010 | [0.016] |  | 0.010 | [0.017]
1 | 0.074 | [0.125] |  | 0.074 | [0.125] |  | 0.074 | [0.125]
8 | 0.400 | [1.000] |  | 0.400 | [1.000] |  | 0.402 | [1.000]
Table 4: Observed maximum convergence rate for the clamped-clamped case, with $n=49$. In brackets are the theoretical rate based on Theorem 3.3. | $p=$
---|---
$K$ | $1$ |  | $2$ |  | $\infty$
$1/8$ | 0.0003 | [0.0033] |  | 0.0003 | [0.00033] |  | 0.0003 | [0.00033]
1 | 0.0020 | [0.0026] |  | 0.0020 | [0.0026] |  | 0.0020 | [0.0026]
8 | 0.0158 | [0.0209] |  | 0.0159 | [0.0209] |  | 0.0161 | [0.0209]
32 | 0.0615 | [0.0836] |  | 0.0619 | [0.0836] |  | 0.0627 | [0.0836]
128 | 0.2223 | [0.3344] |  | 0.2237 | [0.3344] |  | 0.2262 | [0.3344]
Table 5: Observed maximum convergence rate for the clamped-clamped case, with $n=100$. In brackets are the theoretical rate based on Theorem 3.3 | $p=$
---|---
$K$ | $1$ |  | $2$ |  | $\infty$
$1/8$ | 0.0002 | [0.00033] |  | 0.0002 | [0.00033] |  | 0.0003 | [0.00033]
1 | 0.0020 | [0.0026] |  | 0.0020 | [0.0026] |  | 0.0020 | [0.0026]
8 | 0.0157 | [0.0208] |  | 0.0159 | [0.0208] |  | 0.0160 | [0.0208]
32 | 0.0614 | [0.0834] |  | 0.0618 | [0.0834] |  | 0.0625 | [0.0834]
128 | 0.2218 | [0.3336] |  | 0.2232 | [0.3336] |  | 0.2257 | [0.3336]

##  5 Conclusion

The explicit inverse formula for pentadiagonal matrices arising in the fourth-
order nonlinear beam boundary value problem were constructed. The explicit
formula helped computing some norms of their inverse, used to estimate the
convergence of a fixed-point iteration for solving the nonlinear system of
equations. Further research on the convergence upper bounds is necessary to
extend our knowledge of the convergence rate in the fixed point method.

## Acknowledgment

BK and YA wishes to acknowledge the research grant, No AP08052762, from the
Ministry of Education and Science of the Republic of Kazakhstan and the
Nazarbayev University Faculty Development Competitive Research Grant
(NUFDCRG), Grant No 110119FD4502.

## References

  * [1] Chaojie Wang, Hongyi Li, and Di Zhao.  An explicit formula for the inverse of a pentadiagonal toeplitz matrix.  Journal of Computational and Applied Mathematics, 278:12–18, 2015\.
  * [2] EL Allgower.  Exact inverses of certain band matrices.  Numerische Mathematik, 21(4):279–284, 1973.
  * [3] Wayne W Barrett and Philip J Feinsilver.  Inverses of banded matrices.  Linear Algebra and its Applications, 41:111–130, 1981.
  * [4] Lars Rehnqvist.  Inversion of certain symmetric band matrices.  BIT Numerical Mathematics, 12(1):90–98, 1972.
  * [5] Murray Dow.  Explicit inverses of toeplitz and associated matrices.  ANZIAM Journal, 44:185–215, 2002.
  * [6] Wayne W Barrett.  A theorem on inverse of tridiagonal matrices.  Linear Algebra and Its Applications, 27:211–217, 1979.
  * [7] Jiteng Jia, Tomohiro Sogabe, and Moawwad El-Mikkawy.  Inversion of k-tridiagonal matrices with toeplitz structure.  Computers & Mathematics with Applications, 65(1):116–125, 2013\.
  * [8] P Rózsa.  On the inverse of band matrices.  Integral equations and operator theory, 10(1):82–95, 1987.
  * [9] Xi-Le Zhao and Ting-Zhu Huang.  On the inverse of a general pentadiagonal matrix.  Applied Mathematics and Computation, 202(2):639–646, 2008.
  * [10] F Diele and L Lopez.  The use of the factorization of five-diagonal matrices by tridiagonal toeplitz matrices.  Applied mathematics letters, 11(3):61–69, 1998.
  * [11] Xiao-Guang Lv and Ting-Zhu Huang.  A note on inversion of toeplitz matrices.  Applied Mathematics Letters, 20(12):1189–1193, 2007.
  * [12] Bakytzhan Kurmanbek, Yogi Erlangga, and Yerlan Amanbek.  Inverse properties of a class of seven-diagonal (near) toeplitz matrices.  arXiv preprint arXiv:2103.09868, 2021.
  * [13] Mohamed Elouafi.  Explicit inversion of band toeplitz matrices by discrete fourier transform.  Linear and Multilinear Algebra, 66(9):1767–1782, 2018.
  * [14] Yerlan Amanbek, Zhibin Du, Yogi Erlangga, Carlos M. da Fonseca, Bakytzhan Kurmanbek, and António Pereira.  Explicit determinantal formula for a class of banded matrices.  Open Mathematics, 18(1):1227–1229, 2020.
  * [15] Bakytzhan Kurmanbek, Yerlan Amanbek, and Yogi Erlangga.  A proof of andjelić-fonseca conjectures on the determinant of some toeplitz matrices and their generalization.  Linear and Multilinear Algebra, pages 1–8, 2020.
  * [16] Yaroslav Shitov.  The determinants of certain (0, 1) toeplitz matrices.  Linear Algebra and its Applications, 2021.
  * [17] H Grob.  Schwelldruck im belchentunnel.  In Proc. Int. Symp. für Untertagebau, Luzern, pages 99–119, 1972.
  * [18] PA Von Wolffersdorff and S Fritzsche.  Laboratory swell tests on overconsolidated clay and diagenetic solidified clay rocks.  Proc. Geotechnical Measurements and Modelling, Karlsruhe, AA Balkema Pub, pages 407–412, 2003.
  * [19] Piotr Skrzypacz, Daulet Nurakhmetov, and Dongming Wei.  Generalized stiffness and effective mass coefficients for power-law euler–bernoulli beams.  Acta Mechanica Sinica, pages 1–16, 2019.
  * [20] Dongming Wei, Yu Liu, Dichuan Zhang, Match Wai Lun Ko, and Jong R Kim.  Numerical analysis for retaining walls subjected to swelling pressure.  In Proceedings of 2016 International Conference on Architecture, Structure and Civil Engineering (ICASCE’16), London (UK) Mar, pages 26–27, 2016\.