AN ASYMPTOTIC APPROACH TO INVERSE SCATTERING PROBLEMS ON WEAKLY NONLINEAR ELASTIC RODS

Elastic wave propagation in weakly nonlinear elastic rods is considered in the time domain. The method of wave splitting is employed to formulate a standard scattering problem, forming the mathematical basis for both direct and inverse problems. A quasi-linear version of the Wendroff scheme (FDTD) is used to solve the direct problem. To solve the inverse problem, an asymptotic expansion is used for the wave field; this linearizes the order equations, allowing the use of standard numerical techniques. Analysis and numerical results are presented for three model inverse problems: (i) recovery of the nonlinear parameter in the stress-strain relation for a homogeneous elastic rod, (ii) recovery of the cross-sectional area for a homogeneous elastic rod, (iii) recovery of the elastic modulus for an inhomogeneous elastic rod.


Introduction
Wave propagation in nonlinear elastic and viscoelastic materials has been an area of interest for some time.There has been an extensive amount of work done on the mathematical modeling of such materials (e.g., [2,10,17,19,20,23]), as well as on the analysis of the behavior of specific materials (e.g., [11,21,26]).These efforts fall into two classes: the modeling of the stress-strain relation, usually in conjunction with curve fitting to experimental data, and the modeling of wave propagation in nonlinear media, both in the time domain and in the frequency domain.A few application areas include seismic imaging for oil exploration, structural dynamics, and nondestructive evaluation.
Extensive work has been done on inverse problems for linear materials [4,5,7,9,12].However, the literature on inverse problems for nonlinear rods is sparse [8,13,18,22,24], especially in the time domain (although there has been recent activity in nonlinear electromagnetic scattering problems using wave splitting [1,6,16,25]).But the time domain is the natural place to study the propagation of short duration pulses, nonperiodic waves and transients, which are important in many applications.Therefore, it is important to develop analytic and numerical tools to examine the time dependent behavior of nonlinear elastic waves.A program of such study was recently begun [8], and further developed in [13].
It is convenient to begin the study of wave propagation in nonlinear elastic media with the study of one-dimensional rods of finite length d, as is done here.This makes the wave splitting analysis tractable, and still provides insight into meaningful applications.

Goals of the paper
At present, it is not clear how best to solve inverse scattering problems in nonlinear rods.An optimization approach was presented in [8].This paper presents an alternative, a simple asymptotic approach for weakly nonlinear elastic rods.There are two goals here.First, the long range intent is to determine, by numerical experimentation, the conditions for which this approach yields reasonable results for some typical inverse problems, in order to provide insight into how to approach inverse problems on any nonlinear rod.Indeed, the authors are beginning an extensive program of study of inverse problems for nonlinear elastic and viscoelastic wave propagation in two-dimensional media, and work is currently underway to improve the numerical implementation of both the optimization and asymptotic approaches, and to search for other feasible methodologies.Second, the immediate intent is to develop an analytical framework for the asymptotic approach so that higher-order corrections can be incorporated into the analysis presented here.It is reasonable to presume that including higher-order corrections could allow the asymptotic approach to resolve inverse problems involving stronger nonlinearities, significantly enhancing the utility of the approach.
The remainder of the paper is organized as follows.Section 2 presents the model problem under consideration.Section 3 outlines the mathematical formulation of scattering problems using wave splitting.The system of equations used in numerical work appears at the end of the section.Section 4 describes the implementation of the direct algorithm.Section 5 contains a brief description of the inverse problems to be considered.In Section 6, there appears the analysis and numerical results for the case of the homogeneous nonlinear elastic rod; the inverse problem is to recover the value of the nonlinear parameter b 1 .Section 7 discusses analysis and numerical results for the case of the homogeneous nonlinear elastic rod with varying cross section.The inverse problem is to recover the cross-sectional area.Section 8 contains analysis and numerical results for the case of the inhomogeneous elastic rod.The inverse problem is to recover the elastic modulus.Section 9 contains a summary of the techniques and results.

Physical model
Consider an infinite nonlinear rod with a slab of nonlinear elastic material lying in x ∈ [0, d].The material properties (density ρ, elastic modulus E, and cross-sectional area A) are constant outside the slab, and vary continuously for all x.In particular, there are no jump discontinuities at x = 0 or x = d.The stress-strain relation is assumed to have the form For purposes of illustration, a quadratic nonlinearity (F(x, ) = E( +b 1

2
)) is presented, although the analysis below could be generalized to other forms of F(x, ) if desired.An important assumption in this paper is that the material is weakly nonlinear, which means that |b 1 | is small.The value of b 1 can be positive or negative, although for most elastic solids, b 1 is negative.
The equation of motion and the compatibility relation are where v(x, t) is particle velocity.This leads to the system of equations where M includes the effects of spatial inhomogeneity and G ultimately affects the speed of waves traveling through the rod in nonlinear fashion (2.5) This system, along with suitable boundary and initial conditions (presented in the next section), forms the mathematical basis for scattering problems on a nonlinear rod.

Problem formulation
In this section, the method of wave splitting is used to put the system (2.4), into a form suitable for numerical computations.This methodology was developed in [8], and is summarized here for convenience.First, write (2.4) in matrix form as It must be assumed that G > 0. For sufficiently small strains, this is not restrictive, since it is reasonable to expect the stress to increase as the strain increases.The system is simplified by defining and introducing new dependent variables

S. Kim and K. L. Kreider 411
For fixed x, g is strictly increasing with respect to , so there exists an inverse g −1 , interpreted as The system is now represented as with N and wave speed c defined by It must be assumed that 1 + 3b 1 ru 1 > 0 for the wave speed to be nonnegative.
The wave splitting transformation is introduced in order to diagonalize the wave speed matrix.New dependent variables are defined by where This final change of variables yields the quasi-linear equations of motion If the material is spatially inhomogeneous, it is advantageous to introduce travel time coordinates in order to straighten the curve for the leading edge in the space-time plane.This makes it easier to set up a discrete grid for numerical work.The speed of the wave front is denoted c τ (x).
The travel time of the wave front for traversing the rod from x = 0 to The travel time coordinate transformation is then At this point, the dynamic equations can be nondimensionalized, dropping tildes, as with c n = c/c τ .Under this scaling, the wave speed along the leading edge is normalized to 1.The explicit forms of the nondimensional coefficients are (3.12) At x = 0, an incident velocity v(0, t) = f(t) is applied, so that the boundary condition is f(t) = −u + (0, t) + u − (0, t).Casuality implies that all fields vanish before the time of first arrival, so that u ± (x, t) = 0 for x > t.This is the formulation used for the direct and inverse scattering problems considered below.

Summary of assumptions
In the analysis presented above, two assumptions are needed to ensure that the formulation is physically meaningful.First, it is necessary that G > 0 for the wave splitting variables to remain real.This means that 1 + 2b 1 (x, t) > 0, which implies that the strain field remains small in magnitude.The value of can be positive (the rod is in extension) or negative (the rod is in compression).If b 1 > 0, then > −1/2b 1 , which means the rod can be in extension or slight compression.If b 1 < 0, then < −1/2b 1 = 1/2|b 1 |, which means the rod can be in compression or slight extension.Second, it is necessary that (1 + 3b 1 ru 1 ) > 0 for the wave speed to remain nonnegative.The split field The assumption that |b 1 | is small is not necessary here, but does make these inequalities more easily satisfied; the assumption is used in later sections to develop asymptotic algorithms for inverse problems.

The direct problem
A quasi-linear version of the Wendroff scheme [3] is applied to the dynamic equations (3.11), in which the partial derivatives are discretized as For linear equations (c n is constant), the scheme is second order, with error O(∆x 2 + ∆t 2 ), and is unconditionally stable.There is no theory to predict the stability properties of the quasi-linear form, but no numerical difficulties have appeared in the problems under consideration here during extensive testing.This discretization leads to an implicit scheme, the result of which is a coupled system along each time line t j .The unknowns {u ± i,j , i = 1, . . ., j − 1} are obtained one row at a time.Leading edge values u ± jj are computed earlier, as described below.Figure 4.1 shows the grid schematic for an arbitrary row.The dynamic equations are applied at the cell centers, indicated by an ×, providing 2j − 2 equations.The system is closed by including the boundary condition f j = −u + 0j + u − 0j and the directional derivative of u − at the leading edge, obtained from the dynamic equation for u The presence of the nonlinear factors N and c n complicates the analysis; these factors should be evaluated at cell centers (x i+1/2 , t j+1/2 ), presumably by averaging the values of u + and u − at the four surrounding grid points.But this leads to a nonlinear system along time line t = t j because u ± i,j and u ± i+1,j are unknown.Since the material is assumed to be weakly nonlinear, this difficulty can be avoided by averaging the known values of u + and u − at the two lower grid points (i, j − 1) and (i + 1, j − 1).
The quasi-linear direct algorithm is used to create synthetic reflection and transmission data to be used as inputs in the inverse problems presented below.

Three inverse problems
Three inverse problems considered in this paper are the following.
(1) The homogeneous elastic rod.The cross-sectional area A(x) = 1, the constant density ρ, and constant elastic modulus E are known.The goal is to recover the nonlinear parameter b 1 .
(2) The rod with varying cross section.The density, modulus, and nonlinear parameter are known.The goal is to recover the cross-sectional area A(x).
(3) The rod with varying modulus.The density, cross section, and nonlinear parameter are known.The goal is to recover the modulus E(x).
Figure 5.1 contains pictures of the rods and the computational domains for each problem.
In general, we would expect these inverse problems for nonlinear materials to be ill-posed.However, if |b 1 | is sufficiently small, then each of the inverse problems has a unique solution.Sketches of the proofs are presented in the appropriate Sections 6.3, 7.4, and 8.3.Although the details differ from one problem to another, there are two basic themes: the linearized problem has a unique solution if the coefficients are sufficiently smooth and the inputs are sufficiently small, and the higher-order terms in the asymptotic expansion are continuous on the computational domains (Figure 5.1) and so are uniformly bounded.

The homogeneous elastic rod
In this section, the inverse problem for the nonlinear homogeneous elastic rod is formulated, a discussion on uniqueness of solution for the inverse problem is presented, the inverse algorithm is developed, and numerical results are presented and discussed; this material was first described in [13].
Rod geometry for the rod with varying cross-sectional area.(d) Computational domain for the rod with varying crosssectional area and the rod with varying modulus, {(x, t)

Analytic formulation
The nonlinear homogeneous elastic rod is the simplest nonlinear rod.The density ρ and modulus E are constant, so that r(x) can be scaled to 1.Because there are no spatial inhomogeneities, the dynamic equations simplify, considerably, The nonlinear wave speed c n couples the two equations.However, there is no local coupling between u + and u − due to material inhomogeneities, so if no external left-going field (u − ) is applied, then the field is entirely right-going, so u − ≡ 0, and the dynamic equations reduce to a scalar equation in u + as Note that the wave speed c n (u + ) still depends on the field magnitude, and hence the problem is still nonlinear.Also, the boundary condition It is convenient to devise a scattering experiment in which the applied particle velocity at the boundary has the property that f(0) = 0; this implies that u + = 0 along the leading edge, so that the wave front speed c τ = 1 and was used in the numerical experiments presented in Section 6.5.This function was chosen because it offers a smooth, linearly increasing acceleration, which models a linearly increasing force applied to the end of the rod.Because the time interval of interest is [0, 1], there is no issue with the fact that the velocity function is monotonically increasing.
The inverse problem here is to use transmission data, u + (1, t), to recover the nonlinear parameter b 1 .For notational convenience, write u + as u.Straightforward asymptotic expansions of u and c in terms of b 1 , can be applied, resulting in the order equations The important feature in each of these equations is that the wave speed has been linearized; in fact, due to the nondimensionalization, the wave speed is normalized to 1.This means that each field u i propagates along a straight characteristic with slope 1 in space-time coordinates, so that S. Kim and K. L. Kreider 417 the method of characteristics may be used to solve these equations analytically, making the inverse algorithm much more simple.In addition, it should be noted that the cubic term b 1 u 3 was included in the asymptotic expansion, but it was found that this term is so small that it has a negligible effect on the results.Consider (6.4) for u 0 along an arbitrary characteristic t = x + ξ.The equation may be written as which has solution Using ξ = t − x, it is easy to see that partial derivatives of u 0 with respect to x on the curve t = x + ξ may be written as Consider (6.5) for u 1 along the same characteristic t = x + ξ.The equation may be written as which has solution Using ξ = t − x, it is easy to see that partial derivatives of u 1 with respect to x on the curve t = x + ξ may be written as Consider (6.6) for u 2 along the same characteristic t = x + ξ.The equation may be written as which has solution The partial derivatives of u 2 on the curve t = x + ξ may be obtained as above if more terms in the asymptotic expansion of u are desired.
In principle, this process may be repeated to obtain analytic solutions for higher-order terms, although the right side of the order equation becomes increasingly complicated.

Uniqueness of solution theorem
Suppose that the incident velocity v(0, t) = f(t) is smooth and its derivatives are bounded for t ∈ [0, 1].Then for sufficiently small |b 1 |, the inverse problem has a unique solution.

Sketch of proof
Suppose that there are two values of the nonlinear parameter, called b 1 and b 2 , that yield identical transmission data D(t) given the same incident velocity v(0, t) = f(t).Then there are two fields u(x, t) and v(x, t) that satisfy u t + (1 + 3b 1 u) 1/3 u x = 0 and (1, t).The fields may be expanded asymptotically as . From the formulation above, (6.8), (6.11), and (6.14), it is evident that the form of the order functions is independent of the nonlinear parameter, so that u i ≡ v i .It is also evident that the source term in each order equation is continuous, so that each u i = v i is uniformly bounded in the computational domain (Figure 5.1) {(x, t) | 0 ≤ x, x ≤ t ≤ x + 1}.This means that for sufficiently small |b i |, the approximations Because this argument can be carried out to any order of asymptotic expansion, it is clear that b 1 = b 2 , so there is a unique solution to the inverse problem.

Inverse algorithm
The goal of this inverse problem is to recover the nonlinear parameter b 1 .Since there is no reflection in this case, the inverse algorithm is based on transmission data, u + (1, t) = D(t).This data is obtained by experiment or by solving the corresponding quasi-linear direct problem as described above, using the correct value of b 1 as an input.The inverse algorithm is S. Kim and K. L. Kreider 419 based on the least-squares curve fit approximation The transmission data is discretized at time steps t j , and is denoted D j = D(t j ).These data points are added, leading to the sum Because analytic expressions for the u i are given in terms of the known incident field f(t), the only unknown is b 1 , which can be obtained by a standard least squares approach.

Numerical results
A set of numerical experiments was used to gauge the algorithm's overall performance.
Test 1.How much of the transmission data should be used?
Notice that the asymptotic approximation (6.15) to D(t) is quadratic in time, due to the t 2 factor in u 2 , so the least-squares fit will work best in the short time, where D(t) is roughly quadratic.Under the travel time scaling, the wave travels through the material in scaled time t ∈ [0, 1], and transmission data is obtained for this time interval.Let τ denote the fraction of this data that is used by the inverse algorithm.A numerical examination of a variety of cases, with b 1 ranging from −20 to 20, indicates that the best recovery of b 1 occurs when τ < .05(i.e., less than 5% of the data for one travel time is used), although for smaller values of |b 1 |, the value of τ does not appear to be crucial.A reasonable rule of thumb is to let τ = .04for the inverse algorithm to give meaningful results.
Test 2. What is the effect of finite difference discretization error?
The inverse algorithm was run for a variety of b 1 values using τ = .04with a spatial grid of n subintervals.Table 6.The inverse algorithm was run with n = 2048 and τ = .04to determine the range of b 1 values for which the relative error is less than 5% or 10%.To stay within a 10% error bound, the value of b 1 must be in the interval (−12.25,8.75), and to stay within a 5% error bound, the value of b 1 must be in the interval (−5.5, 7.25).
These results clearly indicate that the methodology is useful only for small values of b 1 , and a different technique must be used to analyze strongly nonlinear materials.
Because this problem requires derivatives of the incident field f(t), care must be exercised in specifying a continuous, differentiable incident field.Jumps in f or its derivatives can significantly affect the performance of the algorithm.For example, with the actual value b 1 = −.1, and n = 512 and τ = .40,the algorithm yields the computed value b 1 = −.0838 for the incident field f(t) = −t 2 /2.When this incident field is modified to be equal to .02 for t > .2(so that f is continuous but f has a jump at t = .2),the computed result is b 1 = −.0804, a slight degradation in solution.When a jump is introduced into the incident field (f(t) = t 2 /2 for t < .2 and f(t) = .5for t ≥ .2), the algorithm yields b 1 = 1.831, which is clearly unacceptable.
Overall, the three tests indicate that the recovery of b 1 is reasonable for small values of b 1 , but that this algorithm has a limited usefulness in general.

An inverse problem for the elastic rod with varying cross section
In this section, the inverse problem for the nonlinear homogeneous elastic rod with varying cross-sectional area is formulated, the inverse algorithm is developed, a discussion of uniqueness of solution for the inverse problem is presented, and numerical results are presented and discussed; this was first described in [13].

Analytic formulation
For this problem, the stress-strain relation is again taken to be where |b 1 | < 1 and the scalings ρ = 1, E = 1 are used.The dynamic equations take the form where α(x) = (1/2)A (x)/A(x), and the nonlinear wave speed c n and nonlinear source term F couple the two equations.The leading terms of the asymptotic expansion of u + and u − satisfy As in the previous example, the wave speed has been linearized, but the equations remained coupled.This means that the inverse algorithm for the homogeneous rod cannot be used, because analytic expressions for u ± 0 are no longer available.A different approach must be taken.

Inverse algorithm
Consider the recovery of A(x) from reflection data u − (0, t) = D(t) without knowledge of b 1 .The reflection data is obtained by solving the direct problem, using the correct values of b 1 and A(x) as inputs.Then, the leading order equations (7.3) and (7.4) are solved numerically to obtain the values of A(x i ) at each spatial grid point.The main source of error in this algorithm, which is described below, is that the input data D(t) has not been linearized.Therefore, the method is appropriate only for weakly nonlinear materials, for which the correction terms (u 1 , u 2 ,...) may be safely ignored.
For notational convenience, denote u ± 0 as u ± .The inverse problem is specified as follows: the dynamic equations (7.3) and (7.4) are combined with boundary conditions and leading edge conditions where f(t) = v(0, t) is the known incident field and D(t) is the reflection data, to form a well-posed problem.Note that it is necessary that f(0) = 0 because the leading edge condition is used to recover A(x).For this reason, the initial velocity v(0, t) = f(t) = 1 + t was chosen for this case, as a matter of convenience.
The leading edge conditions are derived from a propagation of discontinuities argument.First, the curve t = x is not a characteristic for u − , so u − cannot admit a discontinuity along that curve.Since u − = 0 for x > t by causality, it must therefore be zero along the leading edge.However, t = x is a characteristic for u + , so there may be a discontinuity in u + along the leading edge.Consider (7.3) for t = x + and t = x − , and subtract the two to get an equation for the jump in u + , denoted Along t = x, this is an ordinary differential equation with [u − ] = 0, which has solution given by (7.6).
The inverse algorithm is a finite difference approach to the method of characteristics.First, discretize the domain by letting ∆x = h, ∆t = 2h, with h = 1/N for some integer N, so that x i = (i − 1)h for i = 1, . . ., N + 1, and t j = 2(j − 1)h for j = 1, . . ., N + 1.Note that t j is measured in wave front time; that is, it is the time after the first arrival of the wave front.

S. Kim and K. L. Kreider 423
The absolute time at spatial location x i is given by x i + t j .The grid points are numbered so that the leading edge corresponds to j = 1, with higher values of j representing later right-moving characteristics, given by t = x + t j .This grid should be consistent with that used in the direct problem, so that the values D(t j ) are easily accessible.
To obtain equations valid at grid location (i, j), apply forward differencing at (i − 1, j) and backward differencing at (i, j) to the u + equation, (7.3), and add, to obtain Then apply forward differencing at (i − 1, j + 1) and backward differencing at (i, j) to the u − equation, (7.4), and add, to obtain Taking values at spatial location i − 1 to be known, (7.9) and (7.10) form a 2 × 2 system in unknowns u + i,j and u − i,j , which can be easily solved to give where (7.12) Equation (7.11) is valid at any interior point in the grid (not on the boundary i = 1 or the leading edge j = 1).
To determine the leading edge conditions at j = 1, use (7.4).Apply backward differencing at (i, 1) to u − equation to obtain where α i = (1/2)A i /A i and A i = (A i − A i−1 )/h.Since u − 0 = 0, then (7.13) becomes All terms are known except A i , therefore, A i can be solved numerically in the equation below using Newton's method.
The inverse algorithm proceeds as follows.
(1) The value of A 1 has been scaled to 1, so begin with the back characteristic t = t 2 − x, denoted by k = 2.There are two grid points along this line, one on the boundary i = 1 and one on the leading edge j = 1, so the boundary condition and leading edge equation are used to obtain A 2 .
(2) Move up to the next back characteristic t = t k − x for k = 3, 4, . . ., N + 1.The boundary condition provides values for u + and u − at i = 1, so use (7.11) for i = 2, 3, . . ., k − 1 to obtain u + and u − at the interior grid points along the characteristic.Then for i = k, the leading edge, (7.15) is used to obtain A k .The algorithm proceeds until k = N + 1, so that the cross-sectional area is obtained at each grid point.Data D(t) is required for t ∈ [0, 2], so that the incident wave can travel through the rod and reflect back from the far end.

Uniqueness of solution theorem
Suppose that α(x) = A (x)/2A(x) is smooth, that |b 1 | is sufficiently small, and that the reflection data D(t) = u − (0, t) is sufficiently small in magnitude.Then the inverse problem, to recover A(x) from the reflection data, has a unique solution.

Sketch of proof
Suppose that there exist two smooth cross-sectional areas A(x) and B(x) that yield the same reflection data The key step is to establish that the linear problem, (7.3) and (7.4), has a unique solution under appropriate conditions.Such a result has been verified for a similar linear system.In [14,15], it is shown that if the reflection and transmission data are sufficiently small, then there is a unique solution to the inverse problem of recovering C(x) and D(x) from that data for u xx − u tt + C(x)u x + D(x)u t = 0.It is assumed that C and D are continuously differentiable.When wave splitting is applied S. Kim and K. L. Kreider 425 to this equation, the result is the system 3) and (7.4) has a similar but simpler form; in particular, it has only one unknown parameter, so transmission data is not required.At any rate, the argument used in [15] can be used successfully here to verify that the inverse problem for the linear system, (7.3) and (7.4), has a unique solution as long as α is smooth and the reflection data D(t) is sufficiently small in magnitude.
Suppose that an initial velocity f(t) is applied to the rod at x = 0, and that cross-sectional areas A(x) and B(x) yield the fields u ± (x, t) and v ± (x, t), respectively, as solutions to (7.2).Expand both fields asymptotically to obtain u ).Then u ± 0 and v ± 0 both satisfy (7.3) and (7.4) with α being replaced by α 0 and β 0 , respectively.Let w ± (x, t) = u ± (x, t)− v ± (x, t), and expand w ± as w ).Then w ± 0 satisfies the linear system Also, w − 0 (0, t) = 0 because both A(x) and B(x) yield the same reflection data, and w + 0 (0, t) = 0 because both u ± and v ± are generated using the same initial velocity.By the uniqueness of the linear problem, w ± 0 (x, t) ≡ 0, so u ± 0 (x, t) = v ± 0 (x, t).Then, since v + 0 + v − 0 ≡ 0, the system for w ± 0 reduces to 0 = ±(α 0 − β 0 )(v + 0 v − 0 ), so α 0 (x) = β 0 (x).A similar argument is used to show that α 1 (x) = β 1 (x).The only difference here is that the source for u ± 1 , v ± 1 includes ∂ x u ± 0 , ∂ x v ± 0 , respectively, so it is necessary to have u ± 0 and v ± 0 sufficiently smooth to be able to claim that u ± 1 and v ± 1 are continuous, which in turn is needed to be able to invoke the uniqueness of the linear problem.Since α(x) and β(x) are smooth, it is true that u ± 0 and v ± 0 are smooth inside the computational domain.
The argument can be extended to the higher-order asymptotic terms; for |b 1 | sufficiently small, the asymptotic approximation of u ± can be made arbitrarily close to the actual value, so that the uniqueness of the linear problem can be extended to the weakly nonlinear case.
It should be noted that the condition that the reflection data be small is sufficient but not necessary for uniqueness-in fact, an explicit example is given in [15].As is often the case, the numerical algorithms perform much better in practice than the theory guarantees.

Numerical results
The inverse algorithm was tested using a variety of cross-sectional area profiles, with n = 128 and b 1 = ±.01,.1,.2.The profiles were chosen with     are not included because the results are similar in nature to those presented.
Using n = 128 grid points in these examples is sufficient.Testing indicates that increasing the number of grid points has very little effect on the results.A typical case is shown in Figure 7.4, where A(x) = 1 − .1sin(2πx)with b 1 = +.1.It should be noted that the inverse algorithm is virtually instantaneous, taking less than 1 second of CPU time on a Pentium III at 600 Mhz, using the g77 compiler on a FORTRAN code in Red Hat Linux.The generation of synthetic data using the data algorithm takes at most a few minutes.
Overall, the results are not as good as we would like, because the inverse algorithm requires very small magnitudes of b 1 to provide reasonable results.However, the inverse algorithm is extremely fast, and does provide good results for weakly nonlinear rods.

An inverse problem for the elastic rod with varying modulus
In this section, the inverse problem for the nonlinear elastic rod with varying modulus of elasticity is formulated, the inverse algorithm is developed, a discussion of uniqueness of solution for the inverse problem is presented, and numerical results are presented and discussed.S. Kim and K. L. Kreider 429

Analytic formulation
For this problem, the stress-strain relation is again taken to be where b 1 1 and the scalings ρ = 1, A = 1 are used.The dynamic equations take the form where, for this case, 3) The nonlinear source term F, as well as the wave speeds, couple the two equations.The leading order asymptotic expansion of the dynamic equations yields It is convenient to apply the travel time coordinate transformation which converts (8.4) into the computational forms Here, H(z) = E(x(z)), H = dH/dz, and v ± (z, s) = u ± (x, t).
Again, the wave speed has been linearized and the equations remain coupled.The dynamic equations may be solved using the method of characteristics.

Inverse algorithm
Consider the recovery of H(z) from reflection data v − (0, s) = D(s) without knowledge of b 1 .The reflection data is obtained by solving the direct problem, using the correct values of b 1 and H(z) as inputs.Then, the leading order equations (8.6) are solved numerically to obtain the values of H(z i ) at each spatial grid point.Then the travel time transformation is used to obtain the value of physical value of x i associated with each computational value z i .As with the varying cross-sectional area, the method is appropriate only for very weakly nonlinear materials.
The inverse problem is specified as follows: the dynamic equations (8.6) are combined with boundary conditions and leading edge conditions where f(s) is the known applied velocity at the left boundary and D(s) is the reflection data, to form a well-posed problem.Equation (8.8) is obtained in the same manner as (7.6).At this point, the inverse algorithm is the same as that for the varying cross-sectional area, with one difference.The algorithm provides discrete values of H i at the points z i .These values must be converted to physical coordinates to obtain the desired E(x i ) values.These values are obtained by discretizing the travel time transformation equation and solving iteratively for x i : (8.10)

Uniqueness of solution theorem
Suppose that H (z)/4H 3/2 (z) is smooth, that |b 1 | is sufficiently small, and that the reflection data D(s) = u − (0, s) is sufficiently small in magnitude.Then the inverse problem, to recover A(x) from the reflection data, has a unique solution.

Inverse scattering on weakly nonlinear rods
The proof follows the format used to show uniqueness for the varying cross-sectional area problem in the previous section.

Numerical results
A number of numerical tests were performed to study this case.The following modulus profiles were studied: The results are shown in Figures 8.1 and 8.2.The main observation to be made here is that the inverse algorithm is much more sensitive to variations in modulus than it is to variations in the cross-sectional area.This makes sense, because the modulus appears in the stress-strain relation and hence plays a more significant role in determining the nature of wave propagation through the rod.However, the results match qualitatively with those for the varying cross-sectional area: small |b 1 | values lead to good reconstructions, larger |b 1 | values lead to worse reconstructions, and if b 1 is large enough, the algorithm fails because the wave speed becomes negative.

Conclusion
In this paper, elastic wave propagation in weakly nonlinear elastic rods is considered in the time domain.The method of wave splitting is employed to formulate a standard scattering problem, forming the mathematical basis for both direct and inverse problems.The focus here is on developing algorithms for solving the inverse problem.An asymptotic approach is used to linearize the dynamic equations.The goal is to determine the conditions for which this approach is reasonable, and to use the insight gained here to develop more robust inverse algorithms.
For the homogeneous nonlinear rod, the asymptotic approach leads to an analytic expression for the transmitted field, which is used in a least squares sense to recover the nonlinear parameter b 1 .For the homogeneous rod with varying cross section and the rod with varying modulus of elasticity, the method of characteristics is used to recover the respective material parameter as a function of depth into the rod.
Numerical results indicate that although the inverse algorithms are extremely fast, the asymptotic approach yields good results only when the nonlinearity in the stress-strain relation is very weak.For even moderate nonlinearities, another approach is needed for solving inverse problems.
Despite the limited usefulness of the algorithms presented here, there are some valuable insights to be gained from this work.The elastic modulus affects wave propagation much more strongly than does the crosssectional area in nonlinear elastic rods, because of its appearance in the stress-strain relation.The results here indicate that it is reasonable to pursue a method that works in a global sense-these algorithms work well near the incident boundary, with results that gradually worsen deeper into the rod, and with occasional sudden failures.An algorithm that does not use this "layer stripping" approach may work much better.
This paper is a first step in determining practical algorithms for solving inverse problems on nonlinear rods in the time domain.Work is currently underway to include higher-order asymptotic terms, to investigate the effects of noisy data on this method, and to consider other approaches to solving such inverse problems.One such approach is the optimization method presented in [8].Although the context here is elasticity, the analysis and numerics work just as well in electromagnetics and acoustics.

Call for Papers
Thinking about nonlinearity in engineering areas, up to the 70s, was focused on intentionally built nonlinear parts in order to improve the operational characteristics of a device or system.Keying, saturation, hysteretic phenomena, and dead zones were added to existing devices increasing their behavior diversity and precision.In this context, an intrinsic nonlinearity was treated just as a linear approximation, around equilibrium points.
Inspired on the rediscovering of the richness of nonlinear and chaotic phenomena, engineers started using analytical tools from "Qualitative Theory of Differential Equations," allowing more precise analysis and synthesis, in order to produce new vital products and services.Bifurcation theory, dynamical systems and chaos started to be part of the mandatory set of tools for design engineers.
This proposed special edition of the Mathematical Problems in Engineering aims to provide a picture of the importance of the bifurcation theory, relating it with nonlinear and chaotic dynamics for natural and engineered systems.Ideas of how this dynamics can be captured through precisely tailored real and numerical experiments and understanding by the combination of specific tools that associate dynamical system theory and geometric tools in a very clever, sophisticated, and at the same time simple and unique analytical environment are the subject of this issue, allowing new methods to design high-precision devices and equipment.
Authors should follow the Mathematical Problems in Engineering manuscript format described at http://www .hindawi.com/journals/mpe/.Prospective authors should submit an electronic copy of their complete manuscript through the journal Manuscript Tracking System at http:// mts.hindawi.com/according to the following timetable:

Figure 4 . 1 .
Figure 4.1.Numerical grid for the direct algorithm.Circles indicate grid points, and the dynamic equations are applied at cell centers, indicated by ×.

Figure 5 . 1 .
Figure 5.1.Rod geometries and computational domains for the three inverse problems under consideration.(a) Rod geometry for the homogeneous rod and the rod with varying modulus.(b) Computational domain for the homogeneous rod, {(x, t) | 0 ≤ x ≤ 1, x ≤ t ≤ x + 1}.(c) Rod geometry for the rod with varying cross-sectional area.(d) Computational domain for the rod with varying crosssectional area and the rod with varying modulus, {(x, t) | 0 ≤ x, x ≤ t ≤ 2 − x}.
1 shows the computed values of b 1 for a variety of grids for actual values of b 1 = ±.01,±.1, ±1, ±10.The results indicate that finer grids do provide better results, but when computation time is taken into account, the marginal cost is too high.The recommended n value is 2048 for smaller |b 1 | values (roughly speaking, b 1 ∈ (−5, 5)), while it may be advantageous to use n = 4096 if |b 1 | is larger.Test 3. What range of b 1 values can be recovered?

Table 6 . 1 . 10 128 1
The effect of discretization error on the recovery of b 1 .In each section, the actual b 1 value is given above, and the values below are the computed values of b 1 using the indicated number n of spatial grid points.n b 1 = −.01 b 1 = −.1 b 1 = −1 b 1 = −= +.01 b 1 = +.1 b 1 = +1 b 1 = +10

( 7 . 17 )
Results for each profile appear in Figures 7.1, 7.2, and 7.3.Some general observations can be made.First, in each case, the reconstructions are nearly perfect for the smallest values of |b 1 |.As |b 1 | increases, the quality of the reconstructions erodes drastically, although in each case there is an interval starting at x = 0 where the reconstruction is accurate; as |b 1 | increases, this interval shrinks.The reconstructions are worse when b 1 is positive than when it is negative.This is due to the expression for the wave speed c = (1 + 3b 1 ru 1 ) 1/3 /r; since u 1 < 0, it is possible for c to become negative when |3b 1 ru 1 | is large enough.When this happens, the formulation breaks down and the algorithm fails.In the figures, this occurs for b 1 = +.2 in each of the three cases, and for b 1 = +.1 in Figure7.2; the profiles are truncated to avoid needless clutter.Tests were also conducted for A(x) = 1 − .5x,A(x) = 2 − arctan(10)/π− arctan(20(x − .5))/π,and A(x) = 1 − .3sin(2πx).The resulting graphs 1 x) A(x) using n = 128 A(x) using n = 256