PARAMETRIC NONLINEAR PDES WITH MULTIPLE SOLUTIONS: A PGD APPROACH

. This paper presents some insights into the determination, using the Proper Generalized Decomposition, of multiple solutions of nonlinear para- metric partial diﬀerential equations. Although the Proper Generalized Decomposition (PGD) is well suited for computing the solution of, possibly nonlinear, parametric problems that vary smoothly with a physical parameter, no work has been achieved for the case of problems that exhibit multiple solutions for some values of a parameter. For two representative cases, we show how an appropriate parametrization, combined to a nonlinear solution procedure can be devised to describe and compute the multiple solutions of a PDE.


1.
Introduction. Solving computational engineering problems often requires the solution of a partial differential equation (PDE) for many values of a set of parameters as for example in design, optimization or control applications. Despite the availability of efficient dedicated algorithms, softwares and the continuing increase of computational resources, the computational cost associated to these many solutions should not be underestimated, especially in the case of real-time applications or for high fidelity models.
Within the framework of the Proper Generalized Decomposition (PGD) the concept of computational vademecum has been recently introduced to circumvent these shortcomings [6]. One precomputes, in an offline stage, the parametric solution of the problem for a given range of the parameters, to be used online in the context of real-time simulation, optimization, inverse analysis and simulation based control. This is possible through the introduction of the desired parameters as extra-coordinates of the problem. In this approach the use of the PGD is critical as the separated representation of the solution allows circumventing the curse of dimensionality that standard mesh-based discretization techniques suffer when addressing the solution of multidimensional models.
Such parametric solutions have been computed and applied to non-trivial engineering cases for many different types of parameters such as: material, process [9], geometrical [3] or design parameters, initial or boundary conditions [10]. Linear models as well as nonlinear ones have been investigated.
Using the PGD, nonlinear parametric PDEs have been addressed only in the case of a unique parametric solution. To solve these parametric models, the parameter is added as an extra-coordinate of the problem and the PGD is combined with an appropriate linearization technique to overcome the nonlinearity. Different linearization techniques have been successfully investigated [1,14]. Beside the classical Picard and Newton linearizations, the LArge Time INcrement method has been efficiently used to solve many computational mechanics problems [12,18]. Another strategy, the Asymptotic Numerical Method (ANM) [8], has been recently used together with the PGD to compute non-linear parametric solutions of thermal models [13]. Specific advantages of using the ANM are its ability to follow a loading path and efficiently handle bifurcations [8,17,4,5]. Despite the fact that efficient strategies for computing bifurcation diagrams make use of model order reduction methods such as the Proper Orthogonal Decomposition or the Reduced Basis technique [16,11], in our knowledge, no work has considered the PGD applied to nonlinear problems with non unique solution for some values of a parameter.
In this paper we discuss, based on two simple examples, the use of the PGD to capture multiple solutions of nonlinear parametric PDEs. In particular we show that it is possible to compute all the multiple solutions simultaneously if an adequate parametrization is introduced. Additional parameters, considered within the PGD framework as extra-coordinates, allow capturing the solution multiplicity. The introduction of these parameters slightly modifies the original problem (by perturbing the original PDE) or the solution procedure (by considering a suitable algorithm for solving the nonlinear problem). More details on both approaches will be given later.
The paper is organized as follows. In section 2, the PGD is briefly revisited. In section 3, we discuss the parametric solutions of nonlinear PDEs with multiple solutions. The strategy is illustrated for two problems: a beam-buckling case and a non-symmetric Fisher-like equation. Section 4 presents conclusions and perspectives of the present research work.
2. PGD in a nutshell. In this section we briefly recall the PGD in the so-called residual minimization form (see [7] for a more complete introduction to the PGD).
Let us consider a linear differential problem with appropriate boundary conditions, whose solution u belongs to an appropriate functional space. We consider a multiparametric case where the solution writes u(x, λ, s), with (x, λ, s) ∈ Ω × I × S, where x is the space coordinate and λ and s are two parametric coordinates. The PGD, applied to this problem, consists in computing the solution approximation u h of u under the following finite sum (also called separated representation): Denoting u h n a solution composed of n terms, the next approximation u h n+1 is such that: is usually computed using a fixed point method, where the minimization of the above expression is carried out in turn on each of the three functions X n+1 (x), L n+1 (λ) and S n+1 (s). The interested reader can refer to the primer [7] for the details on the separated representation constructor.
In this way the PGD algorithm only requires the solution of smaller problems defined on the lower dimensional spaces Ω, I and S, thus avoiding the curse of dimensionality. This enrichment process is repeated until the norm of the discrete residual becomes small enough or convergence is reached according to any other appropriate indicator [2].
3. Nonlinear problems with multiple solutions. Let us consider a parametric problem whose solution u(x, λ) is not unique when the parameter λ is greater than some critical value λ c . In practice, whether in the real world or in simulations, the observation of a specific solution occurs through a perturbation of the original problem which selects a particular solution among all the possible ones. Such perturbation can have different origins and meanings. It could for example represent some material defect or a slightly different initial condition in computational mechanics models or simply a particular initialization in a linearization scheme.
Within the PGD framework, we suggest the introduction of an additional coordinate s which parametrizes a modification of the original problem in such a way that the modified problem no longer has multiple solutions. Indeed each of the original multiple solutions is now the unique solution of the modified problem for a specific value of the additional coordinate s. In what follows, we illustrate this strategy on two simple examples.
3.1. A beam buckling case. We study the compression of a beam of length l clamped at one of its extremities (x = 0) and subject to a compression force F at the other end (x = l). If the force is small enough the problem has a unique solution for which the beam stays aligned in the force direction and experiences only compression. Beyond a critical force F c , this solution becomes unstable and two stable buckled solutions appear as well. Depending on the characteristics of the beam, buckling occurs at a critical load denoted λ c . The theoretical critical load is: where E is the Young modulus, F c is the critical applied force, and I z is the area moment of the beam cross section.
This model problem has been largely studied and several numerical techniques allow the computation of the different bifurcation branches. As explained above, we assume a small perturbation of the original problem which consists in a small deviation of the applied force with respect to the beam axis direction. Let us denote s the angle between the beam and the force direction. The full finite rotation parametric model writes: where θ(x, λ, s) is the angle of the beam axis with respect to the horizontal axis (the undeformed beam direction), and λ = F EIz . The boundary conditions are: θ(x = 0) = 0 and ∂θ ∂x (x = l) = 0, the former describing the clamped condition at x = 0 and the last one the absence of bending moment at x = l. This new model is of course such that, beyond the critical load, the solution θ(x, λ, s) is discontinuous in the s direction around s = 0. Such an expected discontinuity is not an issue for the PGD constructor because the model does not involve derivatives with respect to the added extra-coordinate s.
In order to overcome the numerical pitfalls found close to a bifurcation point, we use the Asymptotic Numerical Method -ANM -to compute the solution where λ is the method loading parameter: 1. The solution, and the parameter λ are expanded into the following power series of the path parameter a: 2. We introduce the above expansions (5,6) into the model equation (4) and we equate the terms with the same power of a. This procedure yields a set of concatenated linear equations for the different orders of the power expansion. 3. Starting from an initial solution θ 0 (x, λ 0 , s), we compute the different unknowns, θ i (x, s) and λ i , of the power expansions. The PGD is then used to compute the solution of each parametric linear problem. The solutions are therefore expressed in the following separated form: For the determination of λ i , we use the following pseudo arc-length criterion: where (·, ·) denotes the scalar product either in L 2 or in R depending on its arguments. 4. After reaching an a priori specified order p in the ANM expansion, we compute its range of validity through: where δ is a small control parameter. 5. If the range of validity is too small, we can perform a continuation step where the whole procedure is restarted. Its starting point being: Remark 1. In the above solution procedure, it should be noted that the use of the PGD for the solution of the parametric linear problem involved at each order of the ANM is critical. Indeed, depending on the case being studied, one could have to solve a multidimensional, multiparametric problem whose solution cannot be addressed by using standard mesh-based discretization techniques.

Remark 2.
Within the same continuation step, the linear operator appearing at each order of the ANM is the same. This feature has made the method very attractive when used with a direct linear solver where the factorization of the discretized system of equations could be recycled from one order to the next. In the case of the PGD solver the assembly of the operator in separated form, which may have a non negligible cost, can therefore be avoided at each order.
Remark 3. The specific form of the ANM solution can easily be recast into the canonical separation of variables used by the PGD. Indeed, as the ANM solution writes: the second expression can be inverted to produce a separated representation of θ(x, λ, s): As θ k (x, s) itself has a separated representation, one can easily rewrite the solution as:  Figure 1 depicts the vertical deflection at the end of the beam with respect to λ and s. As expected, this solution is continuous for λ < λ c and becomes discontinuous with respect to s for s = 0 and λ > λ c . The maximum vertical deflection is again shown in Figure 2 with respect to λ for different values of s: s = 0 • , ±(10 −2 ) • , ±(10 −1 ) • , ±1 • . For very small perturbations (i.e. small s values) the method correctly predicts the deflection of the beam around and past the critical load. The first two stable branches of the bifurcation diagram actually correspond to the two one-sided limits of the solution around s = 0.
To illustrate the performance of the ANM on this parametric problem, we plot in Figure 3 the maximum value of lambda attained after each ANM step. We observe that close to λ c , the size of each step dramatically diminishes.
In order to compute all solutions (u + (x, λ) and u − (x, λ)) simultaneously, we solve Eq. (15) with an additional discrete coordinate s ∈ {−1, 1} which controls the initialization of the iterative procedure from: As the problem now becomes multiparametric we use a PGD solver at each iteration of Eq. (15), where we look for the solution under the following separated form: where Ω x = [0, 1], I λ = [λ min , λ max ] and s ∈ {−1, 1}.
Remark 4. In this second example, the parameter s has no physical meaning. This parameter, considered as an extra-coordinated within the PGD framework, is a purely algorithmic parameter which allows to consider all initial iterates of the linearization procedure.
Remark 5. The solution of Eq. (15) requires the evaluation of the cube of the previous iterate u 3 j . Since the solution is searched in separated form, the evaluation of this term can be time consuming but it does not represent a major difficulty.
3.2.1. Numerical results for the non-symmetric Fisher-like equation. The converged PGD solution u(x, λ, s) of the problem presented in the previous section is obtained with 40 terms in its finite sum representation. Convergence is reached after approximately 150 iterations. Figure 4 depicts u(x, λ, s) for s = 1 and s = −1. For λ < λ c the solution does not depend on s, but for λ > λ c we can clearly observe a positive and a negative part, corresponding to u + and u − respectively. The iterates u j (x = 0.5, λ, s = ±1) are plotted in Figure 5 for different values of j. The dashed line indicates the value of λ c . We see that the convergence is quite fast, except when λ ≈ λ c . In this area, as many as 150 iterations are needed for reaching convergence.

Conclusion.
We have shown in two simple examples that the PGD can be used to compute parametric solutions of nonlinear problems in cases where the combined effect of the parameter value and the nonlinearity can yield multiple solutions. In such cases, it is not sufficient to consider only the parameter as an extra-coordinate of the problem but it is also necessary to introduce an additional coordinate in order to parametrize the different solutions. In order to ensure that the multiple solutions are well parametrized by this additional coordinate we have proposed two choices. First, one can relate the additional coordinate to a physical perturbation of the model that will select a particular solution. In our case the extra-coordinate was the direction of the applied load controlling beam buckling. Second we considered an additional coordinate of discrete nature parametrizing the initialization of the nonlinear solver. In both cases some a priori knowledge of the problem is necessary, but when the appropriate choice is made one can obtain in one shot the different solution branches.