Using Numerical Experiments To Discover Theorems In

Preparing to load PDF file. please wait...

0 of 0
Using Numerical Experiments To Discover Theorems In

Transcript Of Using Numerical Experiments To Discover Theorems In


Website: pp. 495–504

JOSEPH D. FEHRIBACH Department of Mathematical Sciences
Worcester Polytechnic Institute Worcester, MA 01609-2247
Abstract. This work explores the use of numerical experiments in two specific cases: (1) the discovery of two families of exact solutions to the elastic string equations, and approximately periodic solutions that appear to exist near pseudo-solutions formed from these families; (2) the study of the diffusion-reaction-conduction process in an electrolyte wedge (meniscus corner) of a current-producing porous electrode. This latter work establishes the well-posedness of the electrolyte wedge problem and provides asymptotic expansions for the current density and total current produced by such a wedge. The theme of this paper is the use of computing to discover a result that is difficult or impossible to find without a computer, but which once observed, can then be proven mathematically.
1. Introduction. Over the past twenty years, numerical experiments have become an increasingly important part of work in applied mathematics. The current article presents two examples of results which were discovered using computers, then proved using mathematics.
The first example pertains to a classical problem of a vibrating string and the nonlinear system of equations which describe its motion. The next two sections discuss two families of exact solutions for the hyperbolic equations of this system, and also the approximately periodic solutions to these equations which seem to shadow pseudo-solutions formed from the two families. Section 3 presents a calculation of the total energy associated with these pseudo-solutions; this calculation suggests both when and why the approximately periodic solutions may be particularly stable. Specifically, the formation of shocks in solutions appears to be delayed when the pseudo-solutions nearly preserve energy.
The second example comes from the study of porous electrodes, particularly molten carbonate fuel cell (MCFC) cathodes. Among the many issues that affect the performance of these electrodes, one of the more important is the distribution of electrolyte inside the pores. An interesting aspect of this issue occurs at meniscus corners where the three phases (solid, electrolyte and gas) come together along a curve; the electrolyte wedge problem (EWP) is derived from a two-dimensional cross section of this situation. The version considered here is a steady-state diffusionreaction-conduction problem. The fourth section describes both the proof that the problem is well-posed, and also presents formal asymptotic expansions both for the current density and the total current produced by the wedge. Again for this
1991 Mathematics Subject Classification. 35L65,74J30,35J25. Key words and phrases. numerical experiments, elastic string, exact solutions, periodic pseudosolutions, energy, electrolyte wedge problem, maximum principle, matched asymptotics.



problem, numerical experiments have been helpful if not essential in obtaining the results presented here.
There are certainly many other examples of the use of numerical experiments in applied mathematics, but with the exception of a few early hand computations, almost all come from the last half century, and most from the last twenty years. Their density in time seems to be increasing, and they have important implications for all of mathematics. The computational work of Lorenz in the early 1960’s, for example, might have been more readily appreciated had mathematicians of that time been more disposed to consider numerical experiments [8, 9].

2. Elastic String Equations. Consider the nonlinear system which governs the motion of an elastic string of length L with fixed endpoints:



(ξ ξ






= rtt,

−1 < x < 1, t > 0,


BC: r(r−(11,, tt)) == ((L0,,00)),, t > 0, (1)

IC: rrt((xx,, 00)) == fg((xx)),, −1 < x < 1.

Here r(x, t) gives the position in the plane at time t of the element of the string associated with x ∈ [−1, 1] (the reference configuration); ξ := p2 + q2 is the elongation; (p, q) := rx; and T , the tension, is a material property of the string. Let us assume that T satisfies three conditions (for some ξmax):

(a) T (1) = 0,

(b) T ′(ξ) > T (ξ)/ξ > 0, 1 < ξ < ξmax


(c) T ′′(ξ) < 0,

1 < ξ < ξmax

It is often convenient to recast the PDE of (1) as

Ut = F(U)x


where U := (p, q, u, v) := (rx, rt), and F(p, q, u, v) := (u, v, pT (ξ)/ξ, qT (ξ)/ξ). Numerous authors have studied this or similar systems, e.g., Mihailescu & Suli-
ciu [10], Liu [7], Antman [1], Shearer [13, 14, 15] and Rosenau & Rubin [12]. Keller & Ting [6] established that no purely transverse (i.e., u ≡ 0) periodic solutions to (1) exist, but it remains unclear if any periodic solution with nonzero longitudinal motion exist, or how close to periodic an actual solution can be. Shearer [15] solved the Riemann problem for (3) under the assumption that tension satisfies (2). Using this solution and the method of Glimm [5], Fehribach & Shearer [4] ran a variety of numerical experiments and observed surprisingly stable oscillations (approximately periodic solutions) for certain initial data1. In addition, the numerical experiments suggested that these approximately periodic solutions are “close” to two families of special solutions whose existence had not previously been suspected. These observations lead one to the following theorem [4]:

1Interestingly the original reason for running these experiments was to study shocks, not their absence. Although Glimm’s method has only first-order accuracy, it does represent shocks as sharp jumps, and this was why it was chosen.



Theorem 2.1. Equation (3) admits two families of solutions with zero horizontal velocity:
(i) If U = Ub(x, t) = (p(x, t), q(t), 0, v(x)) is a C1 solution of (3), then there are constants p0, q0, v0 and c such that

Ub(x, t) = (p0, q0 + ct, 0, v0 + cx).


(ii) If U = Um(x, t) = (p(x, t), q(x), 0, v(t)) is a C1 solution of (3), then there are constants α, v1 and β > 0 such that

Um(x, t) = (Q(x; α, β), αxQ(x; α, β), 0, v1 + αβt)



T −1 β 1 + (αx)2

Q(x; α, β) :=



1 + (αx)2

Proof. (i) Substituting the ansatz U = Ub(x, t) = (p(x, t), q(t), 0, v(x)) into (3), one immediately obtains that pt = ux and qt = vx, implying the given forms for q and v and that pt = 0. That px = 0 also follows since

T (ξ)

T (ξ) p2

ξ p x = T (ξ) − ξ ξ2 px = ut = 0, (7)

and the bracketed expression is positive because of the tension assumptions (2).

(ii) Again one must substitute the ansatz U = Um(x, t) = (p(x, t), q(x), 0, v(t)) into (3); following the implications again yields the result.

Remark. For the second family, |α| can be thought of as an inverse length scale relative to the reference configuration, while αβ is the acceleration of the string. The solution Ub corresponds to a string swinging about a fixed endpoint (boundary) and contracting, so that the motion of each point on the string is purely vertical. The velocity is increasing in x and constant in time. (cf. Figure 1 (a) ). The solution Um corresponds to an arched middle portion of the string uniformly accelerating in time (cf. Figure 1 (b) ).





Figure 1. The string for the two exact solutions: (a) Corresponding to Ub, the string is swinging about a fixed endpoint and contracting. (b) Corresponding to Um, the string is arched and uniformly accelerating. In this depiction, α < 0.

The special exact solutions of Theorem 2.1 can be combined to form a periodic
pseudo-solution Up of (3) for a quarter period 0 ≤ t ≤ t0 := ξ0/T (ξ0) as follows (cf. Figure 2):

 UL(x, t; a), −1 < x < t/t0 − 1

Up(x, t; a) = UM (x, t; a), t/t0 − 1 < x < 1 − t/t0


 UR(x, t; a),

1 − t/t0 < x < 1



UL(x, t; a) := ξ0(1, a(1 − t/t0), 0, −a(1 + x)/t0) , UM (x, t; a) := (Q(x; a, T (ξ0)), −axQ(x; a, T (ξ0)), 0, −aT (ξ0)t) , (9) UR(x, t; a) := ξ0(1, a(t/t0 − 1), 0, a(x − 1)/t0) .

v q
t = 0 p
v q
t o < t < 2to

v q
0 < t < to p
v q
t = 2t o

p q v
t = to p
q 2t o < t < 3t o








t = 3t o

q 3t o < t < 4t o

q t = 4t o

Figure 2. One complete oscillation for the pseudo-solution Up.
t0 = ξ0/T (ξ0). For the plots which represent time intervals, arrows indicate the motion of the contact points between the two exact solutions. Recall that on each side of a contact point, p, q and v are constant, or depend only on t or x not both, while u ≡ 0 for Up.

Remark. The periodic pseudo-solution Up is defined to satisfy (3) away from the contact points between the two special solutions: xc = ±(1 − t T (ξ0)/ξ0). The velocities and the slope of the string q/p are actually the same for both special solutions at these points, but there is a jump in the elongation and thus the tension. Therefore Up is not itself a weak solution.
2.1. Pseudo-Solution Energy Considerations. Although the pseudo-solution Up defined above is not itself a solution of (3), numerical experiments suggest that for initial data close to Up, no shocks form in the solution of (3) and the solution remains close to Up, at least for several oscillations. These relatively stable oscillatory solutions are referred to as approximately periodic solutions. They appear to have small amounts of longitudinal motion (u = 0) in a neighborhood of each contact point xc between the two exact solutions Ub and Um. This observation is surprising since for continuous initial data far from Up, shocks form quickly, often in less than a single oscillation, and the numerical solutions are not at all



periodic. Thus it is interesting to consider the energy associated with the Up and
the implications of these energy calculations for the elastic string. Let the stored energy function be defined as W (ξ) := 1ξ T (ξ¯) dξ¯, and let the total
energy for the string as a function of time be defined as

E(t) := 1

u2(x, t) + v2(x, t) dx +

W (ξ(x, t)) dx .


2 −1


The first integral in the total energy is the kinetic energy; the second, the potential energy. It is well known that any continuous solution of (1) or (3) must conserve total energy, i.e., dE/dt ≡ 0, while energy must decrease if admissible shocks (discontinuous solutions) occur since such shocks must satisfy an entropy condition (cf. e.g. Shearer [13, p. 451]). Thus for a continuous solution of the full system (1) to in fact be close to the pseudo-solution, Up, the pseudo-solution should at least approximately conserve energy. Hence the question becomes, for reasonable choices of the tension T and the parameters ξ0 and a, is it possible for the approximate solution to nearly conserve energy?
Because of symmetry, one only needs to consider the first quarter of an oscillation: 0 ≤ t ≤ t0 := ξ0/T (ξ0). For convenience, let τ := t/t0; so 0 ≤ τ ≤ 1. Then in terms of τ , the total energy associated with the approximate solution is

E(τ t0) = a2ξ0T (ξ0)τ 2(1 − 2τ /3) + 2τ W (ξ0 1 + a2(1 − τ )2) +



W (T −1(T (ξ0) 1 + a2x2))dx .



Clearly total energy depends on a number of things and is evidently not constant. One can directly compute dE/dτ = t0dE/dt, but the expression is rather complicated, and not particularly enlightening. Using a computer, on the other hand, one can easily show that this energy can be nearly constant for certain parameter values.
At this point, one must consider a specific expression for T (ξ); following the examples in Fehribach & Shearer [4], let T (ξ) := 4 ln(ξ). Using Maple (or similar software), one can quickly plot E in terms of τ . For a = 1 and ξ0 = 2.42, one sees

8.6 19.4
8.59 19.2

8.58 19

8.57 18.8

8.56 18.6

















Figure 3. The total energy E(τ t0) associated with the pseudosolution Up when T (ξ) = 4 ln(ξ), and either a = 1, ξ0 = 2.42 (left) or a = 4, ξ0 = 1.65 (right).

from the plot on the left in Figure 3 that the difference between the minimum and



maximum values of E(τ t0) for τ ∈ [0, 1] is less than 1 % of the minimum value. A similar plot can be obtained when a = 4; in this case, one finds that the optimal choice for ξ0 is approximately 1.65. The right plot in Figure 3 shows that in this case the variation in the total energy associated with the pseudo-solution is again a relatively small percentage of its minimum value. By way of comparison, the numerical experiment shown in Figure 5, Fehribach & Shearer [4] used a = 1 and ξ0 = 2; in that case no discontinuities were observed in the approximately periodic solution through 16384 iterations of Glimm’s method, approximately four complete oscillations of the string (or t ≃ 16t0). This is a relative large number of iterations for a low-order method (Glimm’s method is first order).
The above energy calculations certainly do not establish the existence of periodic solutions to (1). Along with the calculations in Fehribach & Shearer [4], however, they do suggest that there are certain parameter values and tensions T (ξ), and certain initial conditions for which stable continuous oscillatory solutions persist for surprisingly large times, if not indefinitely. None of this, nor Theorem 2.1, would be easily realized without numerical experiments and modern computing.

3. Electrolyte Wedge Problem (EWP). The second case to be considered here comes from the study of mathematical problems relating to the steady-state production of current in porous fuel cell electrodes. These porous electrodes are composed of three phases: the solid electrode phase (Ni or NiO, for example), the electrolyte phase and the gas phase. The latter two phases lie in the void channels (pores) in the solid electrode. In a nutshell, current is produced when the fuel or oxidants diffuse across the gas phase, cross the gas-electrolyte interface, diffuse across the electrolyte phase, and react at the electrolyte-solid interface. The current (electrons or ions) is then conducted out of the solid electrode, through the electrolyte, to a collector plate.
Among the most interesting aspects of this process is what happens in a neighborhood of the triple-contact points (in a cross section perpendicular to a contact curve) where all three phases are in contact. The EWP describes, in terms of the component potentials u and v, the steady-state diffusion, reaction and conduction processes inside this neighborhood. These component potentials are linear combinations of electrochemical potentials for the reacting species; their definition makes possible the study of these processes all on the same scale (Fehribach, Prins-Jansen, et al. [11, 3]). This system is equivalent to the one written in terms of concentrations that is commonly used to model the production of current in such electrodes [2].
The specific problem considered here deals is displayed graphically in Figure 4. This diagram shows a wedge of electrolyte Ωe with contact angle θ0 between a gas phase Ωg and a solid electrode Ωs, a typical situation in MCFC electrodes. The oxidant component potential u is defined in terms of the electrochemical potentials of the oxidants (fuel) that diffuses across Ωg and Ωe, and the current component potential v is defined in terms of the electrochemical potentials of the current carriers in the electrolyte (carbonate ions in the case of a molten carbonate electrode). The only slow (rate-limiting) reaction in this problem occurs at the electrolyte-solid interface ∂Ωes; this slow reaction is represented in the EWP by the ∂Ωes boundary conditions which require that the oxidant potential flux into ∂Ωes equals the current potential flux out of ∂Ωes, and both are proportional to the jump in the component potentials v − u with proportionality coefficient β. The fixed potential difference between the oxidants entering Ωg (u = 0 on ∂Ωg) and the current exiting Ωe (v = v∞ on ∂Ωe) is what drives the production of current in the electrode. The




κ u=0 g
u . n= 0

ε u=0 n κ v=0

u . n= 0 v = v

ε u . n= − κ c v . n= β ( v - u )


Figure 4. Domain Ω := Ωg ∪ Ωe, along with the differential equations and boundary conditions. Note the orientation of the unit vectors at the interfaces. 0 < θ0 < π/2 . The dependent variables u and v are the component potentials associated with the oxidants and the current.

component conductivities κg, ǫ and κc are defined in terms of the oxidant (fuel)

diffusivities and current conductivities (cf. [2] for details). The outer semicircu-

lar boundary is at r = r∞ should be sufficiently distant from the wedge corner


ε β


Written as a system of equations the EWP is


∇ · (κox∇u) = 0

in Ωg ∪ Ωe ,


κc△v = 0

in Ωe ,


∇v · n = 0

on ∂Ωge ,


∇u · n = 0

on ∂Ωsg ,


(e) ε∇u · n = −κc∇v · n = β(v − u) on ∂Ωes ,

(f )



∂Ωg ,

r = r∞

ε β



∇u · n = 0 ; v = v∞


∂Ωe ,

r = r∞

ε β


The notation here is defined as in Figure 4 with the addition that κox := κg in Ωg, while κox := ε in Ωe.
Initially it was suggested that the EWP would be ill-posed in that the current density (which is proportional to ∂nu and ∂nv) would have to be unbounded at the corner. Computational studies using the PDE Toolbox in Matlab with relatively few elements tended to support this view, but larger scale studies suggested otherwise, i.e., that the problem is well-posed. In addition, other computational studies (cf. Figure 5) indicated that, at least to lowest order in ǫ, u is in fact a function of a single similarity variable θ′ defined in the traditional manner relative to a rotated and shifted coordinate system (cf. Figure 6). Given the insights provided by these computational studies, one can prove that the problem is indeed well-posed and



Figure 5. Representative contour plots for v and u in the electrolyte wedge. The plot for v (left) shows the entire electrolyte wedge. Across the entire wedge, v varies almost radially from v∞ (−27020 Joule/mole) on the outer edge (R = 1) to −27017 Joule/mole at the tip of the wedge. Note that near the tip, the contours are not radial with respect to that tip. The plot for u (right) shows only a portion of the wedge near the tip. The contours for u are almost radial lines, but the origin for these lines lies outside the tip of the wedge. The numerical values for u vary from essentially 0 Joule/mole on the upper edge of the wedge to −27020 Joule/mole near the lower right hand corner of the wedge.



Y’ X, X’


x X 0


θ’< 0

Figure 6. Graphical definition of the rotated and shifted coordinate system (the primed coordinate system) in terms of the original coordinate system used in Figure 4. The length scale changes between the lower and upper case coordinate systems by a factor of β/ε; so for example R = rβ/ε.

obtain formal asymptotic expressions for both the current density and the total current.
Theorem 3.1. Let F denote Faraday’s constant, and R∞ := r∞β/ε. The EWP has a unique, bounded solution on the domain Ω. Moreover, the normal derivatives ∂nu and ∂nv (and therefore the current densities) are also bounded along ∂Ωes. Finally formal asymptotic expressions for these normal derivatives can be found in terms of the similarity variable θ′, yielding expressions for both iF , the current



density in the wedge, and its integral, the total current:

iF = β(v − u)/F

= βv∞ θ0 − Tan−1 tan(θ0)

+ O(ǫ, θ′2) ,


F θ0

1 + ε/(θ0βr)


iF dr


εv∞ R∞

θ0 − Tan−1



F θ0

1 + 1/(θ0R∞)

X0 sin(θ0)

R∞ 2


2 ln X0 + 2R∞θ0 + 1


− X0 cos(θ0) θ0 − Tan−1


1 + R∞/(X0 cos(θ0))

The asymptotic matching requires that X0 = cos θ0/θ0.

+ O(ǫ2, θ′2) .

Proof. The proof of the boundedness of the normal derivatives and the uniqueness of the solutions follows from applying the the maximum principle to an equivalent scalar problem where the electrolyte domain for the current component potential v and the reaction interface ∂Ωes are “unfolded” below the domain shown in Figure 4. The unfolding of the domain takes advantage of the decoupling of the equations and boundary conditions for u and v in Ωe. The details of how this equivalent problem is defined are given in [2]. The existence of a solution follows since the equivalent problem is simply a boundary value problem for a generalized Laplace equation with boundary conditions of mixed type. Finally the formal asymptotic expansion for the current density iF follows from a standard (though complicated) application of matched asymptotics (matching outer and inner solutions) to lowest order in both ε and θ′. The asymptotic expansion for the total current follows from direct integration of iF (and was found using the computer algebra system Maple).

Theorem 3.1 has important implications for the material sciences and electrochemistry communities who study the performance of fuel cell electrodes. In electrodes where electrolyte wedges occur, such wedges dominate the production of current. Thus the above asymptotic expressions for current density and total current for an electrolyte wedge can be used along with an estimate for the number of such wedges in a portion of electrode to estimate the overall current density and total current for the electrode. These expressions should also be useful in developing improved homogenized models of the electrodes. Again, these are all results that would not have been found, at least in the author’s opinion, without numerical experiments on modern computers. As the amount of software available continues to increase, along with the power of computer hardware, and as all of the above become increasingly accessible, one can only expect the importance of computing in proving results to continue to increase.
4. Acknowledgment. The author wishes to thank Prof. Schaeffer for long ago introducing him to the power of matched asymptotics.
1. S.S. Antman, The equations for large vibrations of strings, Amer. Math. Monthly 87 (1980), 359–370.
2. J. D. Fehribach, Diffusion-reaction-conduction processes in porous electrodes: the electrolyte wedge problem, European J. Appl. Math. 12 (2001), 77–96.



3. J. D. Fehribach, J. A. Prins-Jansen, K. Hemmes, J. H. W. de Wit, and F. W. Call, On mod-

elling molten carbonate fuel-cell cathodes by electrochemical potentials, J. Appl. Electrochem.

30 (2000), 1015–1021.

4. J. D. Fehribach and M. Shearer, Approximately periodic solutions of the elastic string equa-

tions, Applicable Analysis 32 (1989), 1–14.

5. J. Glimm, Solutions in the large for nonlinear hyperbolic systems of equations, Comm. Pure

Appl. Math. 18 (1965), 95–105.

6. J.B. Keller and L. Ting, Periodic vibrations of systems governed by nonlinear partial differ-

ential equations, Comm. Pure Appl. Math. 19 (1966), 371–420.

7. T.-P. Liu, Development of singularities in the nonlinear waves for quasilinear hyperbolic

partial differential equations, J. Differential Equations 33 (1979), 92–111.

8. E.N. Lorenz, Deterministic non-periodic flow, J. Atmos. Sci. 20 (1963), 130–141.


, The problem of deducing the climate from the governing equations, Tellus 16 (1964),


10. M. Mihailescu and I. Suliciu, Riemann and Gourset step data problems for extensible strings,

J. Math. Anal. Appl. 52 (1975), 20–24.

11. J. A. Prins-Jansen, J. D. Fehribach, K. Hemmes, and J. H. W. de Wit, A three-phase homo-

geneous model for porous electrodes in molten carbonate fuel cells, J. Electrochem. Soc. 143

(1996), 1617–1628.

12. P. Rosenau and M.B. Rubin, Motion of a nonlinear string: Some exact solutions to an old

problem, Phys. Review 31 (1985), 3480–3482.

13. M. Shearer, Elementary wave solutions of the equations describing the motion of an elastic

string, SIAM J. Math. Anal. 16 (1985), 447–459.


, The interaction of transverse waves for the vibrating string, Physical Mathematics

and Nonlinear Partial Differential Equations (J.H. Lightbourne and S.M. Rankin, eds.), Marcel

Dekker, 1985.


, The Riemann problem for the planar motion of an elastic string, J. Differential

Equations 61 (1986), 149–163.

Received Dec. 2002; revised Mar. 2003.
E-mail address: [email protected]