\documentclass[twoside]{article}
\usepackage{amsfonts, amsmath, graphicx}
\pagestyle{myheadings}

\markboth{\hfil An adaptive numerical method for the wave equation
\hfil EJDE/Conf/10} 
{EJDE/Conf/10 \hfil A. S. Ackleh, K. Deng, \& J. Derouen \hfil}

\begin{document}
\setcounter{page}{23}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
Fifth Mississippi State Conference on Differential Equations and
Computational Simulations, \newline
Electronic Journal of Differential Equations,
Conference 10, 2003, pp 23--31. \newline
http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp  ejde.math.swt.edu (login: ftp)}
 \vspace{\bigskipamount} \\
%
  An adaptive numerical method for the wave equation
  with a nonlinear boundary condition
%
\thanks{ {\em Mathematics Subject Classifications:} 35B40, 35L05, 35L20, 65M25, 65N50.
\hfil\break\indent {\em Key words:} Time-adaptive numerical
method, blow-up time, blow-up rate. \hfil\break\indent \copyright
2003 Southwest Texas State University. \hfil\break\indent
Published February 28, 2003. } }

\date{}
\author{Azmy S. Ackleh, Keng Deng, \& Joel Derouen}
\maketitle

\begin{abstract}
 We develop an efficient numerical method for
 studying the existence and non-existence of global solutions to the
 initial-boundary value problem
 \begin{gather*}
 u_{tt}=u_{xx}\quad  0<x<\infty,\;  t>0,\\
 -u_{x}(0,t)=h(u(0,t)) \quad t>0,\\
 u(x,0)=f(x),\quad u_{t}(x,0)=g(x) \quad 0<x<\infty.
 \end{gather*}
 The results by this numerical method corroborate the theory presented in
 \cite{AD}. Furthermore, they suggest that blow-up can occur for more general
 nonlinearities $h(u)$ with weaker conditions on the initial data
 $f$ and $g$.
\end{abstract}

\newtheorem{theorem}{Theorem}[section]
\numberwithin{equation}{section}

\section{Introduction}

In this paper, we consider the initial-boundary value problem
\begin{equation}
\begin{gathered}
u_{tt}=u_{xx}\quad  0<x<\infty, \; t>0,\\
-u_{x}(0,t)=h(u(0,t)) \quad  t>0,\\
u(x,0)=f(x),\quad u_{t}(x,0)=g(x) \quad 0<x<\infty.
\end{gathered} \label{1}
\end{equation}
Here we assume that $h(u)$ is continuous with
$\lim_{u\to\infty}h(u)=\infty$, $g$ is continuous, and
$f$ is continuously differentiable.
To motivate our work for problem (\ref{1}), we point out that this
problem has been recently studied by the authors in \cite{AD}. For
completeness, the main results obtained in that paper are presented as
follows:

\begin{theorem}
There exists at least one mild solution of (\ref{1}) on
$[0,\infty)\times\lbrack0,T_{0})$ for some $T_{0}>0$. Moreover, if $h(u)$
is Lipschitz continuous, then the solution is unique.
\end{theorem}

\begin{theorem}
Suppose that $|h(u)|\leq\rho(|u|)$ with $\rho(r)>0$ continuous, nondecreasing
on $[0,\infty)$, and such that
\[
\int^{\infty}\frac{dr}{\rho(r)}=\infty,
\]
then all mild solutions of (\ref{1}) are global.
\end{theorem}

\begin{theorem}
Suppose that $f(t)+\int_{0}^{t}g(s)ds\geq0\ (\not \equiv0)$ on $[0,\infty)$
and that $h(u)\geq\sigma(|u|)$ with $\sigma(r)>0$ continuous,
nondecreasing on $[0,\infty)$, and such that
\[
\int^{\infty}\frac{dr}{\sigma(r)}<\infty,
\]
then every mild solution of (\ref{1}) blows up in finite time.
\end{theorem}

\begin{theorem}
Suppose that $\int_{0}^{\infty}f(t)dt+\int_{0}^{\infty}\int_{0}^{t}g(s)dsdt>0$
and $h(u)\geq c|u|^{p}$ ($p>1$, $c>0$), then the mild solution of (\ref{1}) blows
up in finite time.
\end{theorem}

In \cite{AD}, we point out that the blow-up occurs on the boundary
$x=0$ only. Moreover, using asymptotic techniques for integral equations
\cite{BH} we establish the following blow-up rates: Letting $T_{b}$ be the
blow-up time,
\begin{itemize}
\item If $h(u)\sim u^{p}$, then $u(0,t)\sim\big(\frac{1}{p-1}\big)
^{\frac{1}{p-1}}(T_{b}-t)^{-\frac{1}{p-1}}$ as $t\to T_{b}$;

\item If $h(u)\sim e^{u}$, then $u(0,t)\sim\log\big( \frac{1}{T_{b}-t}\big)$
as $t\to T_{b}$.
\end{itemize}
The goal of this paper is to develop a numerical method for solving (\ref{1}).
In Section 2 we discuss the numerical approximation while in Section 3, we
present numerical examples. In Section 4, we conclude with some remarks.

\section{Time-Adaptive Method}

We begin this section by integrating (\ref{1}) along characteristics to obtain
the following integral representation of solutions: For $t\leq x$,
\begin{equation}
u(x,t)=\frac{1}{2}[f(x+t)+f(x-t)]+\frac{1}{2}\int_{x-t}^{x+t}g(s)ds,
\label{21}%
\end{equation}
and for $t>x$,
\begin{equation}
\begin{aligned}
u(x,t)=&\frac{1}{2}[f(t+x)+f(t-x)]+\frac{1}{2}\Big[  \int
_{0}^{t+x}g(s)ds+\int_{0}^{t-x}g(s)ds\Big] \\
&+\int_{0}^{t-x}h(u(0,\tau))d\tau.
\end{aligned} \label{22a}
\end{equation}
A solution to the integral equations (\ref{21})-(\ref{22a})
defines a mild solution to the problem (\ref{1}). Furthermore, if
the initial data $f$ and $g$ are smooth and satisfy some
compatibility conditions, then one can argue that a solution of
(\ref{21})-(\ref{22a}) is also a strong solution of (\ref{1}).
Our numerical method will focus on the approximation of (\ref{21}%
)-(\ref{22a}) rather than (\ref{1}). This provides an efficient scheme which
does not require a rather strong regularity assumption on the initial data.

Substituting $x=0$ in (\ref{22a}), we get the Volterra integral equation
\begin{equation}
u(0,t)=f(t)+\int_{0}^{t}g(s)ds+\int_{0}^{t}h(u(0,\tau))d\tau. \label{23}
\end{equation}
Since blow-up occurs only on the boundary $x=0$, a special
attention will be devoted to the development of an approximation
of $u(0,t)$ particularly near the blow-up time $T_{b}$. Once this
is achieved, the approximations of the blow-up time $T_{b}$ and
$u(0,t)$ are used to compute $u(x,t)$ from the equations
(\ref{21})-(\ref{22a}). To this end, differentiating (\ref{23}) we
get the following differential equation for $u(0,t)$:
\[
\frac{du(0,t)}{dt}=\frac{df(t)}{dt}+g(t)+h(u(0,t)).
\]
Let $\Delta t>0$ be sufficiently small. Using Taylor approximation (formally)
we observe that
\[
u(0,t+\Delta t)-u(0,t)=\Delta t\frac{du(0,t)}{dt}+\frac{d^{2}u(0,\xi)}{dt^{2}%
}\Delta t^{2},\quad  \xi\in(t,t+\Delta t).
\]
A key idea in our scheme is to adapt the time step in order to keep the
quantity $|u(0,t+\Delta t)-u(0,t)|\sim|\Delta t\frac{du(0,t)}{dt}|$
bounded by a fixed constant. Since $h(u)\to\infty$ as
$u\to\infty$ and blow-up occurs at $T_{b}$ we see that $
\frac{du(0,t)}{dt}\to\infty$, as $t\to T_{b}$. In particular,
as $t\to T_{b}$ the size of the time step must approach zero if the
magnitude of $\Delta t\frac{du(0,t)}{dt}$ is to remain bounded by a fixed
constant. This forces the numerical approximation not to go beyond the blow-up
time. Making use of this fact we now present a time-adaptive algorithm for
computing $u(0,t)$ and the blow-up time $T_{b}$.

Let $\Delta t_{\min}$ and $\Delta t_{\max}$ be fixed numbers with $0<\Delta
t_{\min}<\Delta t_{\max}<\infty$. Let $u_{0}^{i}$ be the approximation of
$u(0,t_{i})$ with $t_{0}=0$ and $\Delta t_{i}=t_{i}-t_{i-1}\in\left[  \Delta
t_{\min},\Delta t_{\max}\right]$. Denote by
\[
(u_{t})_{0}^{i}=\frac{u_{0}^{i}-u_{0}^{i-1}}{\Delta t_{i}}
\]
the difference approximations of $u_{t}(0,t_{i})$. Guess an
initial time step $\Delta t_{1}$  and fix a scaling factor
$\alpha>1$. Choose constants $d_{l}$ and $d_{u}$ such that
$d_{l}<d_{u}$. The following is a pseudo code for the
time-adaptive algorithm that we have developed:
\\
\texttt{for $i=1,2,\dots$\\
if $\Delta t_{i}|(u_{t})_{0}^{i}|  \leq d_{u}$\\
then\\
 \mbox{\quad} if $i\geq2$\\
 \mbox{\quad} then\\
  \mbox{\quad\quad} if $\Delta t_{i}<\Delta t_{\max}$\\
  \mbox{\quad\quad} then\\
   \mbox{\quad\quad\quad} if $\Delta t_{i}|(u_{t})_{0}^{i}|$
    and $\Delta t_{i-1}|(u_{t})_{0}^{i-1}|\leq d_{l}$ \\
   \mbox{\quad\quad\quad} then \\
    \mbox{\quad\quad\quad\quad}
    $\Delta t_{i+1}=\min(\alpha\,\Delta t_{i},\Delta t_{\max})$\\
   \mbox{\quad\quad\quad} else\\
    \mbox{\quad\quad\quad\quad} $\Delta t_{i+1}=\Delta t_{i}$\\
   \mbox{\quad\quad\quad} end\\
   \mbox{\quad\quad} else\\
   \mbox{\quad\quad\quad} $\Delta t_{i+1}=\Delta t_{i}$\\
  \mbox{\quad\quad} end\\
 \mbox{\quad} else\\
  \mbox{\quad\quad} $\Delta t_{i+1}=\Delta t_{i}$\\
 \mbox{\quad} end\\
 \mbox{\quad} $i=i+1$\\
else\\
 \mbox{\quad} $\Delta t_{i}=\dfrac{\Delta t_{i}}{\alpha}$\\
end \\
done} \medskip

Our adaptive method changes the current time step if one of the
following two cases arises. The first case is that if \textit{\
}$\Delta t_{i}\left|
(u_{t})_{0}^{i}\right|  >d_{u}$ then the approximated quantity $|u_{0}%
^{i+1}-u_{0}^{i}|$ $>d_{u}$. In this case the time step is
decreased by a factor of $1/\alpha$ and the solution is recomputed
at the new time step $(1/\alpha)\Delta t_{i}$. The second case is
that if the current time step $\Delta t_{i}<\Delta t_{\max}$,
$|u_{0}^{i+1}-u_{0}^{i}|\leq d_{l}$ and
$|u_{0}^{i}-u_{0}^{i-1}|\leq d_{l}$, then this indicates that the
time steps used for the last two iterations are very conservative.
Hence, the scheme increases this time step to $\min(\alpha\Delta
t_{i},\Delta t_{\max})$ in order to save computation time. \ It is
easy to see that near the blow-up time, the time step $\Delta
t_{i}$ will decrease until it reaches $\Delta t_{\min}$. When this
happens the computation stops, and the current time is an
approximation of the blow-up time $T_{b}$. We remark that the
accuracy of the approximations of $T_{b}$ depends on the choice of
$\Delta t_{\min}$.

To compute $u_{0}^{i}$ we combine the Runge-Kutta numerical method (see for
example, \cite{FB}) with the above time-adaptive algorithm: Let $u_{0}%
^{0}=f\left(  0\right)  $ and%
\begin{align*}
k_{1} &  =\Delta t_{i+1}y\big(  t_{i},u_{0}^{i}\big)  \\
k_{2} &  =\Delta t_{i+1}y\big(  t_{i}+\frac{\Delta t_{i+1}}{2},u_{0}%
^{i}+\frac{1}{2}k_{1}\big)  \\
k_{3} &  =\Delta t_{i+1}y\big(  t_{i}+\frac{\Delta t_{i+1}}{2},u_{0}%
^{i}+\frac{1}{2}k_{2}\big)  \\
k_{4} &  =\Delta t_{i+1}y\big(  t_{i+1},u_{0}^{i}+k_{3}\big)  ,\,
\end{align*}
where $i=0,1,2,\dots$, and  $\Delta t_{i+1}$ is determined by the
time-adaptive method developed above. Compute $u_{0}^{i+1}$ as follows:
\[
u_{0}^{i+1}=u_{0}^{i}+\frac{1}{6}\left(  k_{1}+2k_{2}+2k_{3}+k_{4}\right).
\]

Now, to approximate the solution of (\ref{21})-(\ref{22a}) we choose $x_{\max
}>0$ and divide the interval $[0,x_{\max}]$ into uniform mesh $x_{j}$ with
$\Delta x=x_{j}-x_{j-1}$, $j=0,1,\dots,m$. Denote by $S^{n}(a,b,I)$ the
Simpson's numerical method for integrating a function $I(t)$ on the interval
$(a,b)$ with $n$ subdivisions, and let $P^{h}(t)$ be the cubic interpolant of
the function $h(u(0,t))$ at the mesh points $t_{i}$. Then we let $u_{j}^{i}$
be the approximation of $u(x_{j},t_{i})$ and compute $u_{j}^{i}$ as follows:
For $t_{i}\leq x_{j}$,
\[
u_{j}^{i}=\frac{1}{2}[f(x_{j}+t_{i})+f(x_{j}-t_{i})]+\frac{1}{2}S^{n}%
(x_{j}-t_{i},x_{j}+t_{i},g),
\]
and for $t_{i}>x_{j}$,
\begin{align*}
u_{j}^{i}   =&\frac{1}{2}[f(t_{i}+x_{j})+f(t_{i}-x_{j})]\\
&  +\frac{1}{2}\left[  S^{n}(0,t_{i}+x_{j},g)+S^{n}(0,t_{i}-x_{j},g)\right]
+S^{n}(0,t_{i}-x_{j},P^{h}).
\end{align*}
In the next section we present numerical results which indicate the accuracy
of such an adaptive numerical scheme in computing both $u(x,t)$ and the
blow-up time $T_{b}$.

\section{Numerical Results}

The numerical method developed in the previous section is now used
to corroborate and complement theoretical results in our earlier
paper \cite{AD}. For the rest of this section let $\Delta
t_{\max}=10^{-3}$, $\Delta t_{\min }=10^{-7}$, $\alpha=2$, $d_u =
1$, $d_l=0.1,$ $n=10,$ $x_{\max}=5$, and $m=200$. In the first
example we present the accuracy of our method. To this end, we
choose $f=0$, $g=1$ and $h(u)=u^{2}$. It is not difficult to show
that $u(0,t)=\tan t,$ and hence blow-up occurs at $t=\pi/2$. In
Figure 1 we show the relative error $\dfrac{|u_{0}^{i}-\tan
t_{i}|}{\Delta t_{i}}$. The computed blow-up time $T_{b}=1.5704$.

\begin{figure}
[ptb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig1.eps}
\caption{The relative error between the computed function $u(0,t)$ and the
exact solution $\tan t$.}
\end{center}
\end{figure}

In the second example we let $f(x)=-(x-2)^{2}+4$, $g(x)=0$ and
$h(u)=u^{3}$. Notice that this choice of initial data does not
satisfy the assumptions of Theorems 1.3-1.4 in Section 1. However,
the numerical results presented in Figures 2-3 indicate that
blow-up occurs for this choice of functions with an approximated
blow-up time $T_{b}=0.5118$.

\begin{figure}[ptb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig2.eps}\caption{The computed function $u(0,t)$ for the data $f(x)=-(x-2)^{2}+4$,
$g(x)=0$ and $h(u)=u^{3}$.}
\end{center}
\end{figure}

\begin{figure} [ptbptb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig3.eps}
\caption{The solution $u(x,t)$ for the data $f(x)=-(x-2)^{2}+4$, $g(x)=0$ and
$h(u)=u^{3}$.}
\end{center}
\end{figure}

In our third numerical experiment we examine whether blow-up
occurs for nonlinearities such as $h(u)=(1+u)[\log(1+u)]^{p}$ with
initial data that do not satisfy the assumptions of Theorem 1.4.
In Figure 4 we present the numerical results of $u(0,t)$ for the
case $p=6$, $f(x)=3e^{-x}\cos(20x)-0.1$ and $g(x)=0$, and in
Figure 5 we display the 3-D graph of the function $u(x,t)$. We
remark that the blow-up time is $T_{b}=0.22296$.

\begin{figure} [ptb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig4.eps}
\caption{The computed function $u(0,t)$ for the data  $f(x)=3e^{-x}%
\cos(20x)-0.1$ and $g(x)=0$ and $h(u)=(1+u)[\log(1+u)]^{6}$.}%
\end{center}
\end{figure}

\begin{figure} [ptb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig5.eps}
\caption{The computed solution $u(x,t)$ for the data
$f(x)=3e^{-x}\cos(20x)-0.1$ and $g(x)=0$ and $h(u)=(1+u)[\log(1+u)]^{6}$.}%
\end{center}
\end{figure}

Using our numerical scheme, we have successfully verified the blow-up rates
given in Section 1 for the functions $e^{u}$ and $u^{p}\,(p>1)$. We now use
this method to examine the blow-up rate for the function $h(u)=(1+u)[\log
(1+u)]^{p}$. Before presenting the numerical results we formally derive such a
rate. Near the blow-up time the values $\frac{df(t)}{dt}$ and
$g(t)$ are negligible when compared to $u(0,t)$, and hence
\[
\frac{du(0,t)}{dt}\sim(1+u(0,t))\left[  \log(1+u(0,t))\right]
^{p}.
\]
Integrating the above we find
\[
\int_{u(0,t)}^{\infty}\frac{du}{(1+u)\left[  \log(1+u)\right]  ^{p}%
}\sim\int_{t}^{T_{b}}dt.
\]
Solving for $u$ we get
\begin{equation}
u(0,t)\sim e^{(\frac{1}{(p-1)(T_{b}-t)})^{\frac{1}{p-1}}}-1.\label{31}%
\end{equation}
In Table 1 we give numerical results that verify such a blow-up rate. For this
computational purpose we use the following equivalent form of (\ref{31})
\[
\frac{1}{p-1}=(T_{b}-t)[\log(1+u(0,t))]^{p-1}.
\]

\begin{table}[ht]
\caption{The blow-up rate for the function $h(u)=(1+u)(\log
(1+u))^{p}$.}
 \begin{center}
\begin{tabular}
[c]{|l|l|l|l|l|}\hline
\multicolumn{1}{|c|}{$p$} & 4 & 6 & 8 & 10\\\hline
Conjectured: $\frac{1}{p-1}$ & $0.3333$ & $0.2$ & $0.1429$ & $0.1111$\\\hline
Approximation & $0.3205$ & $0.1973$ & $0.1411$ & $0.1106$\\\hline
\end{tabular}
\end{center}
\end{table}

\section{Concluding Remarks}

The objective of this paper is to develop a numerical
approximation for studying the existence and non-existence of
global solutions to the wave equation with a nonlinear boundary
condition. Our numerical results indicate that such a scheme is
very accurate and efficient for computing the blow-up time, the
blow-up rate, and the solution. These results also open up several
theoretical questions: 1) How much can the conditions on the
initial data $f$ and $g$ be relaxed for blow-up to occur? 2) Can
one improve Theorem 1.4 for weaker nonlinearties such as
$h(u)=(1+u)[\log(1+u)]^{p}$ ($p>1)$? 3) Can one prove the blow-up
rate  given by (\ref{31}) for such nonlinearities? Our future
research efforts will focus on such questions as well as the
application of time-adaptive methods to a system of wave equations
coupled in the boundary conditions discussed in \cite{AD2}.

Finally, it is worth mentioning that one can also devise a numerical method by
directly approximating the Volterra integral equation (\ref{23}) using a
combination of the time-adaptive method presented here and numerical
quadrature methods for Volterra integral equations \cite{BM}.

\begin{thebibliography}{0} \frenchspacing
\bibitem{AD} A. S. Ackleh and K. Deng, Existence and nonexistence of global
solutions of the wave equation with a nonlinear boundary
condition,\emph{\ Quarterly of Applied Mathematics,} \textbf{59 }(2001), 153-158.

\bibitem{AD2} A. S. Ackleh and K. Deng, Global existence and blow-up for a
system of wave equations coupled in the boundary conditions,
\emph{Dynamics of Continuous, Discrete and Impulsive Systems,}
{\bf 8} (2001), 415-424.

\bibitem{BM} C. T. H. Baker and G. F. Miller, Treatment of Integral Equations by
Numerical Methods, Acedemic Press, New York, 1982.

\bibitem{BH} N. Bleistein and R. A. Handelsman, Asymptotic Expansion of
Integrals, Holt, Rinehart and Winston, New York, 1975.

\bibitem{FB} J. D. Faires and R. L. Burden, Numerical Methods, PWS-Publishing
Company, Boston, 1993.
\end{thebibliography}

\noindent\textsc{Azmy S. Ackleh} (e-mail: ackleh@louisiana.edu)\\
\textsc{Keng Deng} (e-mail: deng@louisiana.edu)\\
\textsc{Joel Derouen} (e-mail: jbd8438@louisiana.edu)\\[2pt]
Department of Mathematics\\
University of Louisiana at Lafayette\\
Lafayette, Louisiana 70504, USA.

\end{document}