Fifth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 10, 2003, pp 23–31. http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

An adaptive numerical method for the wave equation with a nonlinear boundary condition ^{∗}

Azmy S. Ackleh, Keng Deng, & Joel Derouen

Abstract We develop an efficient numerical method for studying the existence and non-existence of global solutions to the initial-boundary value problem

u_{tt} = u_{xx} 0 < x < ∞, t > 0,
−u_{x}(0, t) = h(u(0, t)) t > 0,
u(x, 0) = f(x), u_{t}(x, 0) = g(x) 0 < x < ∞.

The results by this numerical method corroborate the theory presented in [1]. Furthermore, they suggest that blow-up can occur for more general nonlinearities h(u) with weaker conditions on the initial data f and g.

1 Introduction

In this paper, we consider the initial-boundary value problem

\begin{matrix}u_{tt} = u_{xx} 0 < x < ∞, t > 0, \\ −u_{x}(0, t) = h(u(0, t)) t > 0, \\ u(x, 0) = f(x), u_{t}(x, 0) = g(x) 0 < x < ∞. \end{matrix}  (1.1)

Here we assume that h(u) is continuous with lim_{u→∞} h(u) = ∞, g is continuous, and f is continuously differentiable. To motivate our work for problem (1.1), we point out that this problem has been recently studied by the authors in [1]. For completeness, the main results obtained in that paper are presented as follows:

Theorem 1.1 There exists at least one mild solution of (1.1) on [0, ∞)×[0, T_{0}) for some T_{0} > 0. Moreover, if h(u) is Lipschitz continuous, then the solution is unique.

Theorem 1.2 Suppose that |h(u)| ≤ ρ(|u|) with ρ(r) > 0 continuous, nondecreasing on [0, ∞), and such that

∫ ^{∞} \frac{dr}{ρ(r)} = ∞,

then all mild solutions of (1.1) are global.

∗Mathematics Subject Classifications: 35B40, 35L05, 35L20, 65M25, 65N50.

Key words: Time-adaptive numerical method, blow-up time, blow-up rate.

\begin{matrix}c \\ ◯ \end{matrix}2003 Southwest Texas State University.

Published February 28, 2003.

Theorem 1.3 Suppose that f(t) + ∫_{0}^{t} g(s)ds ≥ 0 (≢ 0) on [0, ∞) and that h(u) ≥ σ(|u|) with σ(r) > 0 continuous, nondecreasing on [0, ∞), and such that

∫ ^{∞} \frac{dr}{σ(r)} < ∞,

then every mild solution of (1.1) blows up in finite time.

Theorem 1.4 Suppose that ∫_{0}^{∞} f(t)dt + ∫_{0}^{∞} ∫_{0}^{t} g(s)dsdt > 0 and h(u) ≥ c|u|^{p} (p > 1, c > 0), then the mild solution of (1.1) blows up in finite time.

In [1], we point out that the blow-up occurs on the boundary x = 0 only. Moreover, using asymptotic techniques for integral equations [4] we establish the following blow-up rates: Letting T_{b} be the blow-up time,

• If h(u) ∼ u^{p}, then u(0, t) ∼ (\frac{1}{p−1})^{\frac{1}{p−1}} (T_{b} − t)^{− \frac{1}{p−1}} as t → T_{b};

• If h(u) ∼ e^{u}, then u(0, t) ∼ log (\frac{1}{T_{b}−t}) as t → T_{b}.

The goal of this paper is to develop a numerical method for solving (1.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.

2 Time-Adaptive Method

We begin this section by integrating (1.1) along characteristics to obtain the following integral representation of solutions: For t ≤ x,

u(x, t) = \frac{1}{2}[f(x + t) + f(x − t)] + \frac{1}{2} ∫_{x−t}^{x+t} g(s)ds,  (2.1)

and for t > x,

\begin{aligned}u(x, t) = & \frac{1}{2}[f(t + x) + f(t − x)] + \frac{1}{2}[ ∫_{0}^{t+x} g(s)ds + ∫_{0}^{t−x} g(s)ds] \\ & + ∫_{0}^{t−x} h(u(0, τ))dτ. \end{aligned} (2.2)

A solution to the integral equations (2.1)-(2.2) defines a mild solution to the problem (1.1). Furthermore, if the initial data f and g are smooth and satisfy some compatibility conditions, then one can argue that a solution of (2.1)- (2.2) is also a strong solution of (1.1). Our numerical method will focus on the approximation of (2.1)-(2.2) rather than (1.1). This provides an efficient scheme which does not require a rather strong regularity assumption on the initial data.

Substituting x = 0 in (2.2), we get the Volterra integral equation

u(0, t) = f(t) + ∫_{0}^{t} g(s)ds + ∫_{0}^{t} h(u(0, τ))dτ.  (2.3)

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 (2.1)-(2.2). To this end, differentiating (2.3) 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 Δt > 0 be sufficiently small. Using Taylor approximation (formally) we observe that

u(0, t + Δt) − u(0, t) = Δt\frac{du(0, t)}{dt} + \frac{d^{2}u(0, ξ)}{dt^{2}}Δt^{2}, ξ ∈ (t, t + Δt).

A key idea in our scheme is to adapt the time step in order to keep the quantity |u(0, t+Δt)−u(0, t)| ∼ |Δt\frac{du(0,t)}{dt}| bounded by a fixed constant. Since h(u) → ∞ as u → ∞ and blow-up occurs at T_{b} we see that \frac{du(0,t)}{dt} → ∞, as t → T_{b}. In particular, as t → T_{b} the size of the time step must approach zero if the magnitude of Δ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 Δt_{min} and Δt_{max} be fixed numbers with 0 < Δt_{min} < Δt_{max} < ∞. Let u_{0}^{i} be the approximation of u(0, t_{i}) with t_{0} = 0 and Δt_{i} = t_{i} − t_{i−1} ∈ [Δt_{min}, Δt_{max}]. Denote by

(u_{t})_{0}^{i} = \frac{u_{0}^{i} − u_{0}^{i−1}}{Δt_{i}}

the difference approximations of u_{t}(0, t_{i}). Guess an initial time step Δt_{1} and fix a scaling factor α > 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: for i = 1, 2, … if Δt_{i}|(u_{t})_{0}^{i}| ≤ d_{u}

then

if i ≥ 2

   then

if Δt_{i} < Δt_{max}

     then

if Δt_{i}|(u_{t})_{0}^{i}| and Δt_{i−1}|(u_{t})_{0}^{i−1}| ≤ d_{l}

       then

Δt_{i+1} = min(α Δt_{i}, Δt_{max})

       else

Δt_{i+1} = Δt_{i}

       end
     else

Δt_{i+1} = Δt_{i}

     end
   else

Δt_{i+1} = Δt_{i}

   end

i = i + 1

else

Δt_{i} = \frac{Δt_{i}}{α}

end
done

Our adaptive method changes the current time step if one of the following two cases arises. The first case is that if Δt_{i} |(u_{t})_{0}^{i}| > 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/α and the solution is recomputed at the new time step (1/α)Δt_{i}. The second case is that if the current time step Δt_{i} < Δt_{max}, |u_{0}^{i+1} − u_{0}^{i}| ≤ d_{l} and |u_{0}^{i} − u_{0}^{i−1}| ≤ 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(αΔt_{i}, Δt_{max}) in order to save computation time. It is easy to see that near the blow-up time, the time step Δt_{i} will decrease until it reaches Δ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 Δt_{min}.

To compute u_{0}^{i} we combine the Runge-Kutta numerical method (see for example, [5]) with the above time-adaptive algorithm: Let u_{0}^{0} = f (0) and

k_{1} = Δt_{i+1}y(t_{i}, u_{0}^{i})
k_{2} = Δt_{i+1}y(t_{i} + \frac{Δt_{i+1}}{2} , u_{0}^{i} + \frac{1}{2}k_{1})
k_{3} = Δt_{i+1}y(t_{i} + \frac{Δt_{i+1}}{2} , u_{0}^{i} + \frac{1}{2}k_{2})
k_{4} = Δt_{i+1}y(t_{i+1}, u_{0}^{i} + k_{3}),

where i = 0, 1, 2, …, and Δ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} (k_{1} + 2k_{2} + 2k_{3} + k_{4}) .

Now, to approximate the solution of (2.1)-(2.2) we choose x_{max} > 0 and divide the interval [0, x_{max}] into uniform mesh x_{j} with Δx = x_{j} − x_{j−1}, j = 0, 1, …, 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} ≤ 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),

Figure 1: The relative error between the computed function u(0, t) and the exact solution tan t.

and for t_{i} > x_{j},

u_{j}^{i} =\frac{1}{2}[f(t_{i} + x_{j}) + f(t_{i} − x_{j})]

+ \frac{1}{2} [S^{n}(0, t_{i} + x_{j}, g) + S^{n}(0, t_{i} − x_{j}, g)] + S^{n}(0, t_{i} − x_{j}, P^{h}).

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}.

3 Numerical Results

The numerical method developed in the previous section is now used to corroborate and complement theoretical results in our earlier paper [1]. For the rest of this section let Δt_{max} = 10^{−3}, Δt_{min} = 10^{−7}, α = 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 = π/2. In

Figure 1 we show the relative error \frac{|u_{0}^{i} − tan t_{i}|}{Δt_{i}}. The computed blow-up time

T_{b} = 1.5704.

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 Figure 2: The computed function u(0, t) for the data f(x) = −(x − 2)^{2} + 4, g(x) = 0 and h(u) = u^{3}.

Figures 2-3 indicate that blow-up occurs for this choice of functions with an approximated blow-up time T_{b} = 0.5118.

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.

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 blowup 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} ∼ (1 + u(0, t)) [log(1 + u(0, t))]^{p} .

Integrating the above we find

∫_{u(0,t)}^{∞} \frac{du}{(1 + u) [log(1 + u)]^{p}} ∼ ∫_{t}^{T_{b}} dt.

Solving for u we get

u(0, t) ∼ e^{( \frac{1}{(p−1)(T_{b}−t)} )^{\frac{1}{p−1}}} − 1.  (3.1)

Figure 3: The solution u(x, t) for the data f(x) = −(x − 2)^{2} + 4, g(x) = 0 and h(u) = u^{3}.

Figure 4: 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}.

Figure 5: 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}.

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 (3.1)

\frac{1}{p − 1} = (T_{b} − t)[log(1 + u(0, t))]^{p−1}.

Table 1: The blow-up rate for the function h(u) = (1 + u)(log(1 + u))^{p}.

p  4  6  8  10
Conjectured: \frac{1}{p−1} 0.3333 0.2  0.1429 0.1111
Approximation  0.3205 0.1973 0.1411 0.1106

4 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 (3.1) 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 [2].

Finally, it is worth mentioning that one can also devise a numerical method by directly approximating the Volterra integral equation (2.3) using a combination of the time-adaptive method presented here and numerical quadrature methods for Volterra integral equations [3].

References

[1] A. S. Ackleh and K. Deng, Existence and nonexistence of global solutions of the wave equation with a nonlinear boundary condition, Quarterly of Applied Mathematics, 59 (2001), 153-158.

[2] A. S. Ackleh and K. Deng, Global existence and blow-up for a system of wave equations coupled in the boundary conditions, Dynamics of Continuous, Discrete and Impulsive Systems, 8 (2001), 415-424.

[3] C. T. H. Baker and G. F. Miller, Treatment of Integral Equations by Numerical Methods, Acedemic Press, New York, 1982.

[4] N. Bleistein and R. A. Handelsman, Asymptotic Expansion of Integrals, Holt, Rinehart and Winston, New York, 1975.

[5] J. D. Faires and R. L. Burden, Numerical Methods, PWS-Publishing Company, Boston, 1993.

Azmy S. Ackleh (e-mail: ackleh@louisiana.edu) Keng Deng (e-mail: deng@louisiana.edu) Joel Derouen (e-mail: jbd8438@louisiana.edu) Department of Mathematics University of Louisiana at Lafayette Lafayette, Louisiana 70504, USA.