
We analyze the ground state of a BoseEinstein condensate in the presence of
higherorder interaction (HOI), modeled by a modified GrossPitaevskii equation
(MGPE). In fact, due to the appearance of HOI, the ground state structures
become very rich and complicated. We establish the existence and nonexistence
results under different parameter regimes, and obtain their limiting behaviors
and/or structures with different combinations of HOI and contact interactions.
Both the whole space case and the bounded domain case are considered, where
different structures of ground states are identified.

We propose a new fourthorder compact timesplitting ($S_\text{4c}$) Fourier
pseudospectral method for the Dirac equation by splitting the Dirac equation
into two parts together with using the double commutator between them to
integrate the Dirac equation at each time interval. The method is explicit,
fourthorder in time and spectral order in space. It is unconditional stable
and conserves the total density in the discretized level. It is called a
compact timesplitting method since, at each time step, the number of substeps
in $S_\text{4c}$ is much less than those of the standard fourthorder splitting
method and the fourthorder partitioned RungeKutta splitting method.
Comparison among $S_\text{4c}$ and many other existing timesplitting methods
for the Dirac equation are carried out in terms of accuracy and efficiency as
well as long time behavior. Numerical results demonstrate the advantage in
terms of efficiency and accuracy of the proposed $S_\text{4c}$. Finally we
report the spatial/temporal resolutions of $S_\text{4c}$ for the Dirac equation
in different parameter regimes including the nonrelativistic limit regime, the
semiclassical limit regime, and the simultaneously nonrelativisic and massless
limit regime.

In this paper, we systematically review mathematical models, theories and
numerical methods for ground states and dynamics of spinor BoseEinstein
condensates (BECs) based on the coupled GrossPitaevskii equations (GPEs). We
start with a pseudo spin1/2 BEC system with/without an internal atomic
Josephson junction and spinorbit coupling including (i) existence and
uniqueness as well as nonexistence of ground states under different parameter
regimes, (ii) ground state structures under different limiting parameter
regimes, (iii) dynamical properties, and (iv) efficient and accurate numerical
methods for computing ground states and dynamics. Then we extend these results
to spin1 BEC and spin2 BEC. Finally, extensions to dipolar spinor systems
and/or general spinF (F>=3) BEC are discussed.

We present a regularized finite difference method for the logarithmic
Schr\"odinger equation (LogSE) and establish its error bound. Due to the
blowup of the logarithmic nonlinearity, i.e. $\ln \rho\to \infty$ when
$\rho\rightarrow 0^+$ with $\rho=u^2$ being the density and $u$ being the
complexvalued wave function or order parameter, there are significant
difficulties in designing numerical methods and establishing their error bounds
for the LogSE. In order to suppress the roundoff error and to avoid blowup, a
regularized logarithmic Schr\"odinger equation (RLogSE) is proposed with a
small regularization parameter $0<\varepsilon\ll 1$ and linear convergence is
established between the solutions of RLogSE and LogSE in term of $\varepsilon$.
Then a semiimplicit finite difference method is presented for discretizing the
RLogSE and error estimates are established in terms of the mesh size $h$ and
time step $\tau$ as well as the small regularization parameter $\varepsilon$.
Finally numerical results are reported to confirm our error bounds.

Based on the twodimensional meanfield equations for pancakeshaped dipolar
BoseEinstein condensates in a rotating frame with both attractive and
repulsive dipoledipole interaction (DDI) as well as arbitrary polarization
angle, we study the profiles of the single vortex state and show how the
critical rotational frequency change with the swave contact interaction
strengths, DDI strengths and the polarization angles. In addition, we find
numerically that at the `magic angle' $\vartheta=\arccos(\sqrt{3}/3)$, the
critical rotational frequency is almost independent of the DDI strength. By
numerically solving the dipolar GPE at high rotational speed, we identify
different patterns of vortex lattices which strongly depend on the polarization
direction. As a result, we undergo a study of vortex lattice structures for the
whole regime of polarization direction and find evidence that the vortex
lattice orientation tends to be aligned with the direction of the dipoles.

We study asymptotically and numerically the fundamental gap  the difference
between the first two smallest (and distinct) eigenvalues  of the fractional
Schr\"{o}dinger operator (FSO) and formulate a gap conjecture on the
fundamental gap of the FSO. We begin with an introduction of the FSO on bounded
domains with homogeneous Dirichlet boundary conditions, while the fractional
Laplacian operator defined either via the local fractional Laplacian (i.e. via
the eigenfunctions decomposition of the Laplacian operator) or via the
classical fractional Laplacian (i.e. zero extension of the eigenfunctions
outside the bounded domains and then via the Fourier transform). For the FSO on
bounded domains with either the local fractional Laplacian or the classical
fractional Laplacian, we obtain the fundamental gap of the FSO analytically on
simple geometry without potential and numerically on complicated geometries
and/or with different convex potentials. Based on the asymptotic and extensive
numerical results, a gap conjecture on the fundamental gap of the FSO is
formulated. Surprisingly, for two and higher dimensions, the lower bound of the
fundamental gap depends not only on the diameter of the domain, but also the
diameter of the largest inscribed ball of the domain, which is completely
different from the case of the Schr\"{o}dinger operator. Extensions of these
results for the FSO in the whole space and on bounded domains with periodic
boundary conditions are presented.

We establish uniform error bounds of a finite difference method for the
KleinGordonZakharov system (KGZ) with a dimensionless parameter $\varepsilon
\in (0,1]$, which is inversely proportional to the acoustic speed. In the
subsonic limit regime, i.e. $0<\varepsilon \ll 1$, the solution propagates
highly oscillatory waves in time and/or rapid outgoing initial layers in space
due to the singular perturbation in the Zakharov equation and/or the
incompatibility of the initial data. Specifically, the solution propagates
waves with $O(\varepsilon)$wavelength in time and $O(1)$wavelength in space
as well as outgoing initial layers in space at speed $O(1/\varepsilon)$. This
high oscillation in time and rapid outgoing waves in space of the solution
cause significant burdens in designing numerical methods and establishing error
estimates for KGZ. By adapting an asymptotic consistent formulation, we propose
a uniformly accurate finite difference method and rigorously establish two
independent error bounds at $O(h^2+\tau^2/\varepsilon)$ and
$O(h^2+\tau+\varepsilon)$ with $h$ mesh size and $\tau$ time step. Thus we
obtain a uniform error bound at $O(h^2+\tau)$ for $0<\varepsilon\le 1$. The
main techniques in the analysis include the energy method, cutoff of the
nonlinearity to bound the numerical solution, the integral approximation of the
oscillatory term, and $\varepsilon$dependent error bounds between the
solutions of KGZ and its limiting model when $\varepsilon\to0^+$. Finally,
numerical results are reported to confirm our error bounds.

We present a new approach for predicting stable equilibrium shapes of
crystalline islands on flat substrates, as commonly occur through solidstate
dewetting of thin films. The new theory is a generalization of the widely used
Winterbottom construction i.e., an extension of the Wulff construction for
particles on substrates). This approach is equally applicable to cases where
the crystal surface energy is isotropic, weakly anisotropic, strongly
anisotropic and "cusped". We demonstrate that, unlike in the classical
Winterbottom approach, multiple equilibrium island shapes may be possible when
the surface energy is strongly anisotropic. We analyze these shapes through
perturbation analysis, by calculating the first and second variations of the
total free energy functional with respect to contact locations and island
shape. Based on this analysis, we find the necessary conditions for stable
equilibria and exploit this through a generalization of the Winterbottom
construction to identify all possible stable equilibrium shapes. Finally, we
propose a dynamical evolution method based on surface diffusion mass transport
to determine whether all of the stable equilibrium shapes are dynamically
accessible. Applying this approach, we demonstrate that islands with different
initial shapes may evolve into different stationary shapes and show that these
dynamicallydetermined stationary states correspond to the predicted stable
equilibrium shapes, as obtained from the generalized Winterbottom construction.

We study asymptotically and numerically the fundamental gaps (i.e. the
difference between the first excited state and the ground state) in energy and
chemical potential of the GrossPitaevskii equation (GPE)  nonlinear
Schrodinger equation with cubic nonlinearity  with repulsive interaction
under different trapping potentials including box potential and harmonic
potential. Based on our asymptotic and numerical results, we formulate a gap
conjecture on the fundamental gaps in energy and chemical potential of the GPE
on bounded domains with the homogeneous Dirichlet boundary condition, and in
the whole space with a convex trapping potential growing at least quadratically
in the far field. We then extend these results to the GPE on bounded domains
with either the homogeneous Neumann boundary condition or periodic boundary
condition.

Experiments, theory and atomistic simulations show that finite triple
junction mobility results in nonequilibrium triple junction angles in evolving
polycrystalline systems. These angles have been predicted and verified for
cases where grain boundary migration is steadystate. Yet, steadystate never
occurs during the evolution of polycrystalline microstructures as a result of
changing grain size and topological events (e.g., grain face/edge switching 
"$T_1$" process, or grain disappearance "$T_2$" or "$T_3$" processes). We
examine the nonsteady evolution of the triple junction angle in the vicinity
of topological events and show that large deviations from equilibrium and/or
steadystate angles occur. We analyze the characteristic relaxation time of
triple junction angles $\tau$ by consideration of a pair of topological events,
beginning from steadystate migration. Using numerical results and theoretical
analysis we predict how the triple junction angle varies with time and how
$\tau$ varies with triple junction mobility. We argue that it is precisely
those cases where grain boundaries are moving quickly (e.g., topological
process in nanocrystalline materials), that the classical steadystate
prediction of the finite triple junction mobility triple junction angle is
inapplicable and may only be applied qualitatively.

We study analytically and numerically stability and interaction patterns of
quantized vortex lattices governed by the reduced dynamical law  a system of
ordinary differential equations (ODEs)  in superconductivity. By deriving
several nonautonomous first integrals of the ODEs, we obtain qualitatively
dynamical properties of a cluster of quantized vortices, including global
existence, finite time collision, equilibrium solution and invariant solution
manifolds. For a vortex lattice with 3 vortices, we establish orbital stability
when they have the same winding number and find different collision patterns
when they have different winding numbers. In addition, under several special
initial setups, we can obtain analytical solutions for the nonlinear ODEs.

We propose an efficient and accurate parametric finite element method (PFEM)
for solving sharpinterface continuum models for solidstate dewetting of thin
films with anisotropic surface energies. The governing equations of the
sharpinterface models belong to a new type of highorder (4th or 6thorder)
geometric evolution partial differential equations about open curve/surface
interface tracking problems which include anisotropic surface diffusion flow
and contact line migration. Compared to the traditional methods (e.g.,
markerparticle methods), the proposed PFEM not only has very good accuracy,
but also poses very mild restrictions on the numerical stability, and thus it
has significant advantages for solving this type of open curve evolution
problems with applications in the simulation of solidstate dewetting.
Extensive numerical results are reported to demonstrate the accuracy and high
efficiency of the proposed PFEM.

We present several numerical methods and establish their error estimates for
the discretization of the nonlinear Dirac equation in the nonrelativistic limit
regime, involving a small dimensionless parameter $0<\varepsilon\ll 1$ which is
inversely proportional to the speed of light. In this limit regime, the
solution is highly oscillatory in time, i.e. there are propagating waves with
wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. We
begin with the conservative CrankNicolson finite difference (CNFD) method and
establish rigorously its error estimate which depends explicitly on the mesh
size $h$ and time step $\tau$ as well as the small parameter $0<\varepsilon\le
1$. Based on the error bound, in order to obtain `correct' numerical solutions
in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the CNFD method
requests the $\varepsilon$scalability: $\tau=O(\varepsilon^3)$ and
$h=O(\sqrt{\varepsilon})$. Then we propose and analyze two numerical methods
for the discretization of the nonlinear Dirac equation by using the Fourier
spectral discretization for spatial derivatives combined with the exponential
wave integrator and timesplitting technique for temporal derivatives,
respectively. Rigorous error bounds for the two numerical methods show that
their $\varepsilon$scalability is improved to $\tau=O(\varepsilon^2)$ and
$h=O(1)$ when $0<\varepsilon\ll 1$ compared with the CNFD method. Extensive
numerical results are reported to confirm our error estimates.

We present a uniformly accurate finite difference method and establish
rigorously its uniform error bounds for the Zakharov system (ZS) with a
dimensionless parameter $0<\varepsilon\le 1$, which is inversely proportional
to the speed of sound. In the subsonic limit regime, i.e., $0<\varepsilon\ll
1$, the solution propagates highly oscillatory waves and/or rapid outgoing
initial layers due to the perturbation of the wave operator in ZS and/or the
incompatibility of the initial data which is characterized by two nonnegative
parameters $\alpha$ and $\beta$. Specifically, the solution propagates waves
with $O(\varepsilon)$ and $O(1)$wavelength in time and space, respectively,
and amplitude at $O(\varepsilon^{\min\{2,\alpha,1+\beta\}})$ and
$O(\varepsilon^\alpha)$ for wellprepared ($\alpha\ge1$) and illprepared
($0\le \alpha<1$) initial data, respectively. This high oscillation of the
solution in time brings significant difficulties in designing numerical methods
and establishing their error bounds, especially in the subsonic limit regime. A
uniformly accurate finite difference method is proposed by reformulating ZS
into an asymptotic consistent formulation and adopting an integral
approximation of the oscillatory term. By adapting the energy method and using
the limiting equation via a nonlinear Schr\"{o}dinger equation with an
oscillatory potential, we rigorously establish two independent error bounds and
obtain error bounds at $O(h^2+\tau^{4/3})$ and
$O(h^2+\tau^{1+\frac{\alpha}{2+\alpha}})$ for wellprepared and illprepared
initial data, respectively, which are uniform in both space and time for
$0<\varepsilon\le 1$ and optimal at the second order in space. Numerical
results are reported to demonstrate that our error bounds are sharp.

We analyze rigorously error estimates and compare numerically
spatial/temporal resolution of various numerical methods for the discretization
of the Dirac equation in the nonrelativistic limit regime, involving a small
dimensionless parameter $0<\varepsilon\ll 1$ which is inversely proportional to
the speed of light. In this limit regime, the solution is highly oscillatory in
time, i.e. there are propagating waves with wavelength $O(\varepsilon^2)$ and
$O(1)$ in time and space, respectively. We begin with several frequently used
finite difference time domain (FDTD) methods and obtain rigorously their error
estimates in the nonrelativistic limit regime by paying particular attention to
how error bounds depend explicitly on mesh size $h$ and time step $\tau$ as
well as the small parameter $\varepsilon$. Based on the error bounds, in order
to obtain `correct' numerical solutions in the nonrelativistic limit regime,
i.e. $0<\varepsilon\ll 1$, the FDTD methods share the same
$\varepsilon$scalability on time step: $\tau=O(\varepsilon^3)$. Then we
propose and analyze two numerical methods for the discretization of the Dirac
equation by using the Fourier spectral discretization for spatial derivatives
combined with the exponential wave integrator and timesplitting technique for
temporal derivatives, respectively. Rigorous error bounds for the two numerical
methods show that their $\varepsilon$scalability on time step is improved to
$\tau=O(\varepsilon^2)$ when $0<\varepsilon\ll 1$. Extensive numerical results
are reported to support our error estimates.

We propose and rigourously analyze a multiscale time integrator Fourier
pseudospectral (MTIFP) method for the Dirac equation with a dimensionless
parameter $\varepsilon\in(0,1]$ which is inversely proportional to the speed of
light. In the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the
solution exhibits highly oscillatory propagating waves with wavelength
$O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. Due to the rapid
temporal oscillation, it is quite challenging in designing and analyzing
numerical methods with uniform error bounds in $\varepsilon\in(0,1]$. We
present the MTIFP method based on properly adopting a multiscale decomposition
of the solution of the Dirac equation and applying the exponential wave
integrator with appropriate numerical quadratures. By a careful study of the
error propagation and using the energy method, we establish two independent
error estimates via two different mathematical approaches as
$h^{m_0}+\frac{\tau^2}{\varepsilon^2}$ and $h^{m_0}+\tau^2+\varepsilon^2$,
where $h$ is the mesh size, $\tau$ is the time step and $m_0$ depends on the
regularity of the solution. These two error bounds immediately imply that the
MTIFP method converges uniformly and optimally in space with exponential
convergence rate if the solution is smooth, and uniformly in time with linear
convergence rate at $O(\tau)$ for all $\varepsilon\in(0,1]$ and optimally with
quadratic convergence rate at $O(\tau^2)$ in the regimes when either
$\varepsilon=O(1)$ or $0<\varepsilon\lesssim \tau$. Numerical results are
reported to demonstrate that our error estimates are optimal and sharp.
Finally, the MTIFP method is applied to study numerically the convergence
rates of the solution of the Dirac equation to those of its limiting models
when $\varepsilon\to0^+$.

We propose a sharpinterface continuum model based on a thermodynamic
variational approach to investigate the strong anisotropic effect on
solidstate dewetting including contact line dynamics. For sufficiently strong
surface energy anisotropy, we show that multiple equilibrium shapes may appear
that can not be described by the widely employed Winterbottom construction,
i.e., the modified Wulff construction for an island on a substrate. We repair
the Winterbottom construction to include multiple equilibrium shapes and employ
our evolution model to demonstrate that all such shapes are dynamically
accessible.

We derive rigorously one and twodimensional meanfield equations for cigar
and pancakeshaped BoseEinstein condensates (BEC) with higher order
interactions (HOI). We show how the higher order interaction modifies the
contact interaction of the strongly confined particles. Surprisingly, we find
that the usual Gaussian profile assumption for the strongly confining direction
is inappropriate for the cigarshaped BEC case, and a ThomasFermi type profile
should be adopted instead. Based on the derived mean field equations, the
ThomasFermi densities are analyzed in presence of the contact interaction and
HOI. For both box and harmonic traps in one, two and three dimensions, we
identify the analytical ThomasFermi densities, which depend on the competition
between the contact interaction and the HOI.

We introduce and analyze a novel meanfield model for polariton condensates
with velocity dependence of the effective polariton mass due the photon and
exciton components. The effective mass depends on the inplane wave vector k,
which at the inflection point of the lower polariton energy branch becomes
infinite and above negative. The polariton condensate modes of the new
meanfield theory are now sensitive to mass variations and for certain points
of the energy dispersion the polariton condensate mode represents fractional
quantum mechanics. The impact of the generalized kinetic energy term is
elucidated by numerical studies in 1D and 2D showing significant differences
for large velocities. Analytical expressions for plane wave solutions as well
as a linear waves analysis show the significance of this new model.

A multiscale time integrator Fourier pseudospectral (MTIFP) method is
proposed and analyzed for solving the KleinGordonSchr\"{o}dinger (KGS)
equations in the nonrelativistic limit regime with a dimensionless parameter
$0<\varepsilon\le1$ which is inversely proportional to the speed of light. In
fact, the solution to the KGS equations propagates waves with wavelength at
$O(\varepsilon^2)$ and $O(1)$ in time and space, respectively, when
$0<\varepsilon\ll 1$, which brings significantly numerical burdens in practical
computation. The MTIFP method is designed by adapting a multiscale
decomposition by frequency to the solution at each time step and applying the
Fourier pseudospectral discretization and exponential wave integrators for
spatial and temporal derivatives, respectively. We rigorously establish two
independent error bounds for the MTIFP at $O(\tau^2/\varepsilon^2+h^{m_0})$
and $O(\varepsilon^2+h^{m_0})$ for $\varepsilon\in(0,1]$ with $\tau$ time step
size, $h$ mesh size and $m_0\ge 4$ an integer depending on the regularity of
the solution, which imply that the MTIFP converges uniformly and optimally in
space with exponential convergence rate if the solution is smooth, and
uniformly in time with linear convergence rate at $O(\tau)$ for
$\varepsilon\in(0,1]$ and optimally with quadratic convergence rate at
$O(\tau^2)$ in the regime when either $\varepsilon=O(1)$ or $0<\varepsilon\le
\tau$. Thus the meshing strategy requirement (or $\varepsilon$scalability) of
the MTIFP is $\tau=O(1)$ and $h=O(1)$ for $0<\varepsilon\ll 1$, which is
significantly better than classical methods. Numerical results demonstrate that
our error bounds are optimal and sharp. Finally, the MTIFP method is applied
to study numerically convergence rates of the KGS equations to the limiting
models in the nonrelativistic limit regime.

In this paper, we propose a regularized Newton method for computing ground
states of BoseEinstein condensates (BECs), which can be formulated as an
energy minimization problem with a spherical constraint. The energy functional
and constraint are discretized by either the finite difference, or sine or
Fourier pseudospectral discretization schemes and thus the original infinite
dimensional nonconvex minimization problem is approximated by a finite
dimensional constrained nonconvex minimization problem. Then an initial
solution is first constructed by using a feasible gradient type method, which
is an explicit scheme and maintains the spherical constraint automatically. To
accelerate the convergence of the gradient type method, we approximate the
energy functional by its secondorder Taylor expansion with a regularized term
at each Newton iteration and adopt a cascadic multigrid technique for selecting
initial data. It leads to a standard trustregion subproblem and we solve it
again by the feasible gradient type method. The convergence of the regularized
Newton method is established by adjusting the regularization parameter as the
standard trustregion strategy. Extensive numerical experiments on challenging
examples, including a BEC in three dimensions with an optical lattice potential
and rotating BECs in two dimensions with rapid rotation and strongly repulsive
interaction, show that our method is efficient, accurate and robust.

In this paper, we propose efficient and accurate numerical methods for
computing the ground state and dynamics of the dipolar BoseEinstein
condensates utilising a newly developed dipoledipole interaction (DDI) solver
that is implemented with the nonuniform fast Fourier transform (NUFFT)
algorithm. We begin with the threedimensional (3D) GrossPitaevskii equation
(GPE) with a DDI term and present the corresponding twodimensional (2D) model
under a strongly anisotropic confining potential. Different from existing
methods, the NUFFT based DDI solver removes the singularity by adopting the
spherical/polar coordinates in Fourier space in 3D/2D, respectively, thus it
can achieve spectral accuracy in space and simultaneously maintain high
efficiency by making full use of FFT and NUFFT whenever it is necessary and/or
needed. Then, we incorporate this solver into existing successful methods for
computing the ground state and dynamics of GPE with a DDI for dipolar BEC.
Extensive numerical comparisons with existing methods are carried out for
computing the DDI, ground states and dynamics of the dipolar BEC. Numerical
results show that our new methods outperform existing methods in terms of both
accuracy and efficiency.

We study dimension reduction for the threedimensional GrossPitaevskii
equation with a longrange and anisotropic dipoledipole interaction modeling
dipolar BoseEinstein condensation in a strong interaction regime. The cases of
disk shaped condensates (confinement from dimension three to dimension two) and
cigar shaped condensates (confinement to dimension one) are analyzed. In both
cases, the analysis combines averaging tools and semiclassical techniques.
Asymptotic models are derived, with rates of convergence in terms of two small
dimensionless parameters characterizing the strength of the confinement and the
strength of the interaction between atoms.

We present efficient and accurate numerical methods for computing the ground
state and dynamics of the nonlinear Schr\"odinger equation (NLSE) with nonlocal
interactions based on a fast and accurate evaluation of the longrange
interactions via the nonuniform fast Fourier transform (NUFFT). We begin with a
review of the fast and accurate NUFFT based method in \cite{JGB} for nonlocal
interactions where the singularity of the Fourier symbol of the interaction
kernel at the origin can be canceled by switching to spherical or polar
coordinates. We then extend the method to compute other nonlocal interactions
whose Fourier symbols have stronger singularity at the origin that cannot be
canceled by the coordinate transform. Many of these interactions do not decay
at infinity in the physical space, which adds another layer of complexity since
it is more difficult to impose the correct artificial boundary conditions for
the truncated bounded computational domain. The performance of our method
against other existing methods is illustrated numerically, with particular
attention on the effect of the size of the computational domain in the physical
space. Finally, to study the ground state and dynamics of the NLSE, we propose
efficient and accurate numerical methods by combining the NUFFT method for
potential evaluation with the normalized gradient flow using backward Euler
Fourier pseudospectral discretization and timesplitting Fourier pseudospectral
method, respectively. Extensive numerical comparisons are carried out between
these methods and other existing methods for computing the ground state and
dynamics of the NLSE with various nonlocal interactions. Numerical results show
that our scheme performs much better than those existing methods in terms of
both accuracy and efficiency.

We propose a sharp interface model for simulating solidstate dewetting where
the surface energy is (weakly) anisotropic. The morphology evolution of thin
films is governed by surface diffusion and contact line migration. The
mathematical model is based on an energy variational approach. Anisotropic
surface energies lead to multiple solutions of the contact angle equation at
contact points. Introduction of a finite contact point mobility is both
physically based and leads to robust, unambiguous determination of the contact
angles. We implement the mathematical model in an explicit finite difference
scheme with cubic spline interpolation for evolving marker points. Following
validation of the mathematical and numerical approaches, we simulate the
evolution of thin film islands, semiinfinite films, and films with holes as a
function of film dimensions, Young's angle $\theta_i$, anisotropy strength and
crystal symmetry, and film crystal orientation relative to the substrate
normal. We find that the contact point retraction rate can be well described by
a powerlaw, $l \sim t^n$. Our results demonstrate that the exponent $n$ is not
universal  it is sensitive to the Young's angle $\theta_i$ (and insensitive
to anisotropy). In addition to classical wetting (where holes in a film heal)
and dewetting (where holes in a film grow), we observe cases where a hole
through the film heals but leave a finite size hole/bubble between the
continuous film and substrate or where the hole heals leaving a continuous film
that is not bonded to the substrate. Surface energy anisotropy (i) increases
the instability that leads to island breakup into multiple islands, (ii)
enhances hole healing, and (iii) leads to finite island size even under some
conditions where the Young's angle $\theta_i$ suggests that the film wets the
substrate.