Dang Trinh Nguyen
email: dang-trinh.nguyen@grenoble-inp.fr

Quoc Bao Duong
email: quoc-bao.duong@grenoble-inp.fr

Eric Zamai
email: eric.zamai@grenoble-inp.fr

Muhammad Kashif
email: kashif.shahzad@grenoble-inp.fr

-Nancy -France "Toward circular Economy" Diagnosis in complex system with multiple failure sources

Keywords: Fault diagnosis, complex manufacturing system, Bayesian network

This paper proposes an approach to accurately localize the origin of product quality drifts, in a flexible manufacturing system (FMS). The failure propagation mechanism in a production process is proposed based on the relationships between failure sources to explain the failure propagation following production flow. The logical diagnosis model is used to reduce the search space of suspected equipment in the production flow; however, it does not help in accurately localizing the faulty equipment. In the proposed approach, we model this reduced search space as a Bayesian network that uses historical data to compute conditional probabilities for each suspected equipment. This approach helps in making accurate decisions on localizing the cause for product quality drifts as either one of the equipment in production flow or product itself.

INTRODUCTION

One of main challenges in manufacturing domain is to improve and optimize the production process quality and equipment effectiveness during production operations. To improve product quality and reduce associated costs, fault isolation, detection and diagnosis techniques have been developed. For the fault diagnosis, many methods are proposed through a diagnosis model such as [START_REF] Sampath | Active diagnosis of discrete-event systems, Automatic Control[END_REF][START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF] to localize more quickly and accurately the root causes of a detected failure. Theses methods are appropriate for diagnosis of failures that come from equipment and products. However, these methods do not explain the impact of recipes to product quality. In fact, due to the diversity of customer demands, the production recipes frequently change, and so they may trigger the impacts on product quality.

In addition, the production system comprises of hundreds of equipment, monitored by thousands of sensors. Generally, it is difficult to identify a diagnosis model and corresponding variables. Thus, artificial intelligence methods are well suited techniques to analyze the large amount of production information and describe the behavior of system components as presented in [START_REF] Weber | Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas[END_REF]. These approaches can be performed without understanding the underlying structure of a production system [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF]. Among the probabilistic approaches, the Bayesian network (BN) approach is widely used to identify a graphical structure model that describes relationships between variables in a production system. So, its conditional probabilities will be calculated to provide the risk priorities and support corrective maintenance decisions. Due to the complexity of present-day manufacturing systems, identification of this graphical structure is complex to be performed by a maintenance engineer [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF]. The system elements such as products, recipes, equipment, maintenance schedules and human factors are frequently changed due to the introduction of new production technologies and maintenance management strategies. A change in one element may trigger effect to others. Consequently, we need to retrain the learning phase of the BN model to adapt to newly emerging situations in the production system. The time and work-load for computation are very large, so this result in poor maintenance and additional costs.

Consequently, in this paper, we propose: (i) a failure propagation mechanism to explain the relationships between different failure sources and its propagation in the production process, (ii) a diagnosis methodology that uses the Logical Diagnosis (LD) model [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF] to reduce the search space for faulty equipment from a given production flow and optimizes the learning phase for the subsequent BN. Thus, joint and conditional probabilities of all faulty candidates are computed to support corrective maintenance decisions.

The characteristics of case study and the diagnosis objective are introduced in the next section. Section 3 presents a review on fault diagnosis methods. In Section 4, we propose the failure model that enables to explain impact and propagation of failure in the product flow. The LD model that results in suspected equipment against a failure drift is explained in Section 5. The proposed diagnosis model is described in Section 6. The conclusion and future perspectives are discussed in Section 7.

CASE STUDY AND DIAGNOSIS OB-JECTIVE

Case study

Manufacturing processes evolve to adapt itself to increasing demand diversity. Besides, frequent changes in customer demands, it leads to the changes in corresponding products. This is reason that flexible manufacturing systems (FMS) are widely used in complex and highly competitive manufacturing domains such as (micro-)electronic and automotive industry. The FMS is a complex automated manufacturing system that consists of several production workshops, connected by an appropriate transport system. These production workshops and transport systems are controlled by a control and automation system. It is generally characterized by multiple products, production lines, recipes, and human factors. Generally, a complex system may have many production processes driven by a control law to perform demands of the control system. The structure of a production process describes product type, product lines and corresponding equipment in the operating part of the controlled system. Any change in one of elements (such as product, equipment, recipe and human) through control law leads to changes of the existing processes. Consequently, our research in this paper is based on the context of a general complex system under dynamic production environment within the following characteristics: multiple products, production lines, and recipes as shown in Figure 1. In addition, the failure can come from different sources such equipment, products and recipes.

Diagnosis objective

In complex engineering applications, systems can be composed of many components and subsystems, and the ways in which these elements interact affects the way failures propagate within the subsystems and across subsystem boundaries [START_REF] Hine | A multi-leve/multilayer architecture for intelligent shopfloor control[END_REF]. In this paper, we consider that failure causes may fall among products, equipment and recipes. For monitoring the execution of system components, the hierarchical and modular controls are often used as presented in [START_REF] Hine | A multi-leve/multilayer architecture for intelligent shopfloor control[END_REF]. In this context, an automated manufacturing system is organised by Computer Integrated Manufacturing (CIM) architecture that contains: controlled system, product flow and control system, as shown in Figure 2.

Figure 2: CIM architecture of an automated production system

The controlled system consists of actuators and sensors. The sensors monitor the executions of actuators and product flows. Therefore, actuators and sensors are controlled by a local control module, and the set of these elements are called Functional Chain (FC) [START_REF] Henry | Controlled system model adapted for control law synthesis[END_REF]. The FCs receive demands from the coordination level of control system and executes these demands on product flow. When a FC cannot correctly execute a demand, it implies that a fault is produced. In fact, there are many positions of actuators or products which are not observed by the sensors. Hence, when a failure occurs, it is not observed, it may propagate form one FC to another through production lines; and so it called failure propagation. Therefore, this failure propagation not only causes the potential failure on products, but also may have consequences on other system components; and so, the new faults may be created. Consequently, when a fault is detected by the metrology, its root causes may come from one or a part of elements of system. In summary, the diagnosis objective is to precisely and quickly locate the possible origins of failure to support the maintenance operator to save recovery time (for return to a normal status) of the production system.

LITERATURE REVIEW OF DIAGNO-SIS METHODS

Diagnosis techniques based on failure propagation to analyze component dependencies are propagation graph (Abdelwahed & Karsai, 2006) and temporal chronic [START_REF] Strasser | Diagnostic alarm sequence maturation in timed failure propagation graphs[END_REF]. These methods are based on historical production data to locate components that are possible origins of a detected fault. However, these approaches do not analyse the behaviors of the functional chains (FC), so they cannot explain the consequences between elements of a FC and between different FCs. Moreover, these approaches cannot reduce the size of the model. Due to the cyclical operation of the control system, a large amount of information from the production system provokes the problem of combinatory explosion. Consequently, we are especially interested in the Logical Diagnosis model proposed in [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF]. In this model, a diagnosis function is proposed to characterize the historical information of a controlled system to search the suspected potential fault origins in real time. Hence, this model provides a set of possible origins under the form of a directed graph, and its size is reduced by the exploitation of the controlled system observations. This is appropriate to diagnose the faults of production equipment and products. However, this method can not explain the impact of recipe to product quality and its correlations to other system elements. In fact, any change of a production recipe may trigger a product quality drift. In addition, this model does not show suspected level of each candidate of possible origins in the diagnosis result. It is difficult to decide the maintenance order in a complex system. Consequently, this approach would be extended to optimize corrective maintenance activities.

Concerning the risk priority of potential failure candidates, the probabilistic approaches are widely used such as Neural network [START_REF] Khomfoi | Fault Diagnostic System for a Multilevel Inverter Using a Neural Network[END_REF] and Bayesian network [START_REF] Weber | Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas[END_REF]. These approaches enable to calculate the probability values from the large database and associated variables from a production system. These probability values allow evaluating the suspected level (high or low) in order to support the decision for a maintenance strategy.

In particular, the BN models have the advantages that fit to be applied in manufacturing industry as explained in [START_REF] Weber | Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas[END_REF]. The methods based on BN are introduced in [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF] to make a diagnosis in a multiple variables system. The Confidence level of feedback information is proposed in [START_REF] Duong | Confidence estimation of feedback information for logicdiagnosis[END_REF] to provide the probability value that shows the confidence of faulty execution report from equipment. When the database is available, these approaches are practical tools for the corrective maintenance. However, they must be extended due to the following problems. The structure of confi-dence level model for information feedback (CLFI), as presented in [START_REF] Duong | Confidence estimation of feedback information for logicdiagnosis[END_REF], is static with seven parameters. In the context of flexible manufacturing system with characteristics such as multiple products, production lines, recipes and human factors, the production situation often change as we presented in section 2. Hence, the set of parameters that can have an impact on the confidence level of equipment is dynamic. In addition, the BN model must be updated based on the information of newly production situations. In practice, when the information of database is available, learning approaches are often used for modelling BN as presented in [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]). In the learning approaches, a graphical structure and probabilistic rules are estimated from observed data. Many studies in [START_REF] Bouaziz | Dependability of complex semiconductor systems: Learning Bayesian networks for decision support[END_REF][START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] show that these learning approaches are still complex in identifying variables as they depend strongly on expert opinions. The learning workload for computation is still large [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. This loses too much time and is not appropriate for real time structural identification as it depends on the exploitation of databases. In fact, production environments are increasingly stressed by strong competition. It shows that the time for locating the root causes of failures and process recovery (return the process to a normal status) is very important. These challenges promote the researches to apply BN model for real time fault diagnosis and corrective maintenance optimization.

This paper proposes a diagnosis model that enables real time localization of the possible failure and root causes that come from equipment, products and recipes, thus dynamically computing conditional probabilities between a failure and its possible causes. This is based on a logical diagnosis model and a BN model. Thus, the diagnosis mechanism is constructed based on four main steps: (i) LD model provides a set of possible fault origins. The relationships of members in this set are used to construct the graphical structure, (ii) we use this structure for the BN model. The idea is to simplify the variables identification during learning phase in the BN, (iii) the historical information of production system is used to estimate probabilistic rules in the learning phase of the BN model, (iv) the conditional probabilities of nodes in structural model are computed. These probabilistic values allow to evaluate the risk priority for each possible fault origin. However, the BN model structure depends strongly on cause-consequence relationships between the members in the set of possible origins, while the LD model explains only the propagation of equipment failure and corresponding information (products and recipes), but does not explain their cause-consequence relationships. Consequently, we analyze these relationships and so the failure propagation mechanism in a production process is proposed in the next section.

FAILURE PROPAGATION MECHA-NISM

Relationships between equipment, recipes and product quality

The diagnosis objective in this paper is to accurately and quickly localize the possible failure origins which may come from equipment, products and recipes. In fact, the relationships between failure sources have impacts on final product quality. Thus, these causeconsequence relationships are analyzed in this section. The product quality depends on equipment state as presented in [START_REF] Sampath | Active diagnosis of discrete-event systems, Automatic Control[END_REF]. However, this impact will change when the recipe elements are taken into account. Consequently, we present the definition of recipe and analyze its impacts on equipment and product quality next. The concept of recipe is widely used in the control of Flexible Manufacturing Systems. The representation of the production process is based on the concept of controlled operation with basic and control recipes that are developed through rigorous process R&D used in the automation of batch plants by all major suppliers of process control systems as presented in [START_REF] Mergen | Batch Control Systems Standards[END_REF]. The basic recipe is independent of equipment and this describes the production process for a given technology. It is designed to manufacture the product and its quality depends on the basic recipe. Equipment executes demands of basic recipe through the control recipe. The transition to control recipe consists of the basic recipe, the production schedule, the multi-purpose plant description and executable control recipe [START_REF] Genrich | Verification of recipe-based control procedures by means of predicate/transition nets[END_REF] Generally, the control recipes are created based on the characteristics of system and are tested before applied for the real system [START_REF] Mergen | Batch Control Systems Standards[END_REF]. Hence, the control recipes do not trigger a fault on the equipment. However, different equipment execute the same control recipe and give different product quality due to its different physical characteristics. Moreover, an equipment performance may be good for one recipe, but produces bad quality with other recipes. Con-sequently, there exists the cause-consequence relationships between the equipment, recipe and product quality as shown in Figure 4. It means that product quality drifts may come from one or a part of physical equipment failures and control recipes. In a production process, the failure may propagate from equipment to another through the production flow. Hence, these cause-consequence relationships also propagate following the production flow. This problem is presented in the next section.

Failure propagation mechanism

The failure occurring on equipment can come from physical failure of itself and/or the input products. Therefore, the failure that occurred on equipment and corresponding control recipe may propagate to the out products. In summary, there is causeconsequence relationship between product quality and the set of equipment and recipe. These relation continues to propagate on the product flow, and so they lead to a failure propagation in a production process as shown in Figure 5. The Logical Diagnosis model proposed in [START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF] is one part in the treatment of failure propagation through products in a complex production process. In this model, a diagnosis function is proposed to characterize the historical information from a controlled system in locating the possible origins of a detected fault. This diagnosis function is presented in three main points: First, necessary information is collected from coordination level for fault diagnosis through Operation Models as presented in [START_REF] Henry | Controlled system model adapted for control law synthesis[END_REF]. An Operation Model contains information of Functional Chains (equipment, sensors and local control modular), (Pre-)conditions, (Pre-)constraints and effects of operations as shown in Figure 6. The model for diagnosis describes a graphical structure of a production process that consists of system components (as illustrated in Figure 7) and its relationships following product flows (the arrows → in Figure 7). Second, a mechanism to reduce this model is developed by exploiting controlled system observations. Following an operation in this model, if the information provided by these two elements is not coherent, the controlled system sends to the coordination level a faulty execution report. If the coordination level does not receive a faulty execution report from these Functional Chains, it can conclude that the corresponding element is reliable. These reliable elements (the black nodes in Figure 7) will be removed from the model, while the suspected elements (the white nodes in Figure 7) will be retained.

Finally, a mechanism is defined, based on the failure propagation approach (dash arrows ⇢ in Figure 7), it allows us to search the possible origins and the consequences against a fault when the model receives a faulty execution report from the Functional Chain as described in Figure 7. Figure 7: Structural model For fault diagnosis with equipment failure, products and control recipes, we consider that an Operation O i contains information of an Equipment E i , a corresponding control recipe R i which is a Condition for Operation O i and a product P i associated with Equipment E i which is the final effect of Operation O i . Consequently, we obtain a set of possible faulty equipment {E 1 , ...E i , ...E n }, a set of control recipes {R 1 , ...R i , ...R n } and a set of products {P 1 , ...P i , ...P n }. For instance, a failure propagates from Operation O i to Operation O i+1 . The LD model allows to localize the suspected operation and corresponding equipment, product and recipe. However, it does not describe their relationships inside of the operation. For this reason, we have proposed the failure propagation mechanism in a complex process to explain their cause-consequence as presented in section 4. Hence, the directed graph, that describes the failure propagation from Operation O i to Operation O i+1 within the cause-consequence relationships inside these operations, is represented in Figure 8. Moreover, to evaluate the risk priority of each candidate, the information given by directed graph and sets of equipment, products and recipes will be used for structural identification of the BN model. Therefore the conditional probabilities are computed by the diagnosis model which is presented in the next section.

DIAGNOSIS MODEL

Model description

The proposed model comprises of logical diagnosis model and BN model as shown in Figure 9. In this model, we use the results given by the logical diagnosis model introduced in Section 5 for dynamic structure identification of a BN. After that, the probability values are computed by the BN model to support the decision-making for corrective maintenance. The model for diagnosis is based on the control system and dynamic reconfiguration as proposed in [START_REF] Henry | Controlled system model adapted for control law synthesis[END_REF][START_REF] Deschamps | Diagnosis for control system reconfiguration[END_REF]. The reactive loop is characterized by collaboration of several supervision, monitoring and control (SM&C) functions such as detection, diagnosis, prognosis, decision and automatic control (Zamai et al., 1997). The coordination level of control system has capacity to manage a set of Functional Chains and corresponding information. It demands and memorizes the information from the production database. It also provides the information for diagnosis.

Our methodology consists of: (i) searching possible root causes against a detected fault in the past evolution of the operating part within controlled system, (ii) computing the probability values that show the suspect levels against candidates. The diagnosis model execution is illustrated as following principles:

This model is generated by the coordination level. Thus, it sends commands and receives reports in real time for and from all system components. This coordination level also provides the necessary information for the logical diagnosis.

Once a failure is detected by a metrology, a set of possible origins and its correlations are defined by the logical diagnosis model. This set of possible origins and corresponding information are sent to BN model. Thus, a graphical structure of failure mode is determined to support the structural identification in the learning phase of BN model.

After the establishment of graphical structure for BN model, the conditional probabilities associated with all nodes of network are computed based on historical information in the production database. All computed results are stored in the production database to support decision-making for corrective maintenance.

The execution of diagnosis model

When a fault is detected by a metrology, the diagnosis model manages its execution and will demand the diagnosis results from Logical Diagnosis model and the BN model. The Logical Diagnosis model provides a reduced model of the set of possible faulty equipment {E 1 , ...E i , ...E n }, also the set of corresponding Control recipes {R 1 , ...R i , ...R n } and the set of products {P 1 , ...P i , ...P n } as presented in Section 5. These are sent to BN through coordination level, and used to construct the BN model. The conditional probabilities are computed based on the principal theories as presented in [START_REF] Manfredotti | Modeling and inference with relational dynamic Bayesian networks[END_REF]. The BN model is established and performed as follows:

First, a graphical structure of BN model is transformed from the sets of possible faulty equipment, control recipes and products that are given by logical diagnosis model. Each member in these sets is considered as a node in the BN model. The equipment E i ′ ; (i ′ = 1⋯n, i ′ ≠ i) is the parent of E i if it is in front of the node E i and has direct logical relationship with E i following a product flow in the directed graph. The product P i ′ corresponding with E i ′ is put between E i ′ and E i , while the control recipe against E i ′ is considered a parent of P i ′ .

Consequently, we obtain a graphical structure of the BN model. Each node in this structure may be a parent of child-nodes and may be a child of other parent-nodes. For instance, the parents of a detected fault node are E i and R i . Therefore, E i has a set of parents that are P 1 i-1 , ...P k i-1 , ...O N i-1 , with k j = 1⋯N . Besides, the parents of each node P k i-1 are E k i-1 and R k i-1 . In this case, the BN model has a hierarchical structure as shown in Figure 10. , with k j = 1⋯n j . Hence, the learning phase of BN model is performed to calculate the probabilities P (N j P roduct) for each node N j following the product, and next the conditional probabilities between the child-nodes and its parents based on the historical information of production database. We consider that each member of the set of fault origins has two states {1, 0}. Thus, the conditional probabilities over each node N j with (j = 1⋯n) and its parents are defined as a matrix P N j N 1 j , ...N kj j , ...N nj j

as the next equation:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ p N j = 1 N 1 j = 1, ...N kj j = 1, ...N nj j = 1 ⋮ p N j = 1 N 1 j = 0, ...N kj j = 0, ...N nj j = 0 p N j = 0 N 1 j = 1, ...N kj j = 1, ...N j nj = 1 ⋮ p N j = 0 N 1 j = 0, ...N kj j = 0, ...N nj j = 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (1)
In this paper, in order to compute the conditional probabilities, when we need to extract the distribution over some subset of variables or a single variable, we need to marginalize or sum out the variables other than the variables of interest as explained in [START_REF] Manfredotti | Modeling and inference with relational dynamic Bayesian networks[END_REF]. The marginalization rule for any sets of variables X and Y is given by :

P (X) = y P (X, y); y ∈ Y (2)
The distribution over X can be obtained by summing out all the other variables from any joint distribution containing X. We can use the conditional probabilities instead of joint probabilities to compute the probabilities over X as shown in Equation ( 3 The above sections show that the set of possible origins is dynamically determined when a fault is detected. Therefore, when the data is available, the corresponding conditional probabilities are calculated by the diagnosis model. This set of possible failure origins and these probability values are sent and stored in the production database by the coordination level.

Discussion

The advantages of proposed model are: first to locate the possible fault origin sets in real time, second to reduce the space of this set by the evaluation of suspect levels, and finally less workload for structure identification of the BN model. In the proposed model, the set of possible root causes is significantly reduced by the Logical Diagnosis model. Next, this is used to simplify the structural identification of BN model. Hence, the BN model receives a graphical structure with only elements as possible origins for a detected fault. While other elements not related to the detected fault are removed. Consequently, it does not need to compute all probabilities of all elements in the production system.

In addition, the probability values computed by BN model help to continue reducing the set of possible origins through risk priority to save the recovery time of a production system. This also implies that the combination of a deterministic approach (logical diagnosis) and a probabilistic approach (Bayesian network) help us to locate more accurately and quickly equipment as the root cause of a detected fault. Consequently, this proposed model is feasible to apply for fault diagnosis in complex automated manufacturing systems with large production information.

CONCLUSION

In this paper, the failure propagation mechanism is proposed to describe the impact of recipe on equipment and product quality, and to explain the failure propagation mechanism in a production process. Hence, the proposed diagnosis model dynamically generates the structure of the BN and the associated probabilities. We used a Logical Diagnosis model to significantly reduce the search space for suspected equipment in the given production flow. This reduced set of possible origins which is demonstrated as the directed graph provides the cause-consequent relation to simplify the failure model identification in the learning phase of BN. In addition, the associated probabilities are computed by the BN model to evaluate the suspect level of each member in the set of possible fault origins. Consequently, The proposed model is appropriate in dynamically locating the root causes, in less time and less workload of the computation of conditional probability values in the context of complex manufacturing system that is characterized by multiple products, production lines and recipes. Thus, the diagnosis results support decision-making for corrective maintenance activities.

In future work, we are interested in applying the proposed diagnosis method to production systems to evaluate performances of the proposed method. Therefore, the corresponding algorithms need to be developed in order to adapt to real time manufacturing systems. In addition, we will improve the fault diagnosis in a general manufacturing system with multiple failure sources.

Figure 1 :

 1 Figure 1: A Flexible Manufacturing System

  as illustrated in Figure 3.

Figure 3 :

 3 Figure 3: The process to Create a Control Recipe

Figure 4 :

 4 Figure 4: The relationship between different failure sources

Figure 5 :

 5 Figure 5: Failure propagation in a production process

Figure

  Figure 6: Operation Model

  Indeed, this logical diagnosis model provides a reduced structural model and a set of suspected operations that have logical relationships with faulty execution. These suspected operations are considered as possible origins {O 1 , ...O i , ...O n } against detected fault. The reduced model describes the logical links between possible fault origins; thus, these logical links are considered as the cause-consequence relationships according to the failure propagation. The form of reduced model is a directed graph in which the nodes represent suspected operations, while the arcs represent the paths of failure propagation through product flows.

Figure 8 :

 8 Figure 8: Possible fault origins inside operations

Figure 9 :

 9 Figure 9: Diagnosis model

Figure 10 :

 10 Figure 10: Graphical structure of BN model

  model computes conditional probabilities P (N j f ault, P roduct) over nodes N j with (j = 1⋯n) given by detected fault according to the Product as illustrated in the next equation:P (N j f ault, P roduct) = n i=1 Ni P (N j ).P (N i N j ).P (f ault N 1 , ..., N i , ..., N n , P roduct) (4) with i = 1⋯n, i ≠ j. They are used to support the maintenance decision.


Martin J Gander
email: martin.gander@unige.ch

Laurence Halpern
email: halpern@math.univ-paris13.fr

AND V Éronique

A NEW ALGORITHM BASED ON FACTORIZATION FOR HETEROGENEOUS DOMAIN DECOMPOSITION

Keywords: Heterogeneous domain decomposition, viscous problems with inviscid approximations, transmission conditions, factorization algorithm AMS subject classifications. 65M55, 65M15

Often computational models are too expensive to be solved in the entire domain of simulation, and a cheaper model would suffice away from the main zone of interest. We present for the concrete example of an evolution problem of advection reaction diffusion type a heterogeneous domain decomposition algorithm which allows us to recover a solution that is very close to the solution of the fully viscous problem, but solves only an inviscid problem in parts of the domain. Our new algorithm is based on the factorization of the underlying differential operator, and we therefore call it factorization algorithm. We give a detailed error analysis, and show that we can obtain approximations in the viscous region which are much closer to the viscous solution in the entire domain of simulation than approximations obtained by other heterogeneous domain decomposition algorithms from the literature.

1. Introduction. The coupling of different types of partial differential equations is an active field of research, since the need for such coupling arises in various applications. A first main area is the simulation of complex objects, composed of different materials, which are naturally modeled by different equations; fluid-structure interaction is a typical example, and many techniques have been developed for this type of coupling problems, see for example the book [START_REF] Morand | Fluid Structure Interaction[END_REF], or the review on the immersed boundary method [START_REF] Mittal | Immersed boundary methods[END_REF], and [START_REF] Deparis | Heterogeneous domain decomposition methods for fluid-structure interaction problems[END_REF] for domain decomposition coupling techniques. A very important area of application is the simulation of the cardiovascular system [START_REF] Formaggia | Cardiovascular Mathematics: Modeling and simulation of the circulatory system[END_REF]. A second main area is when homogeneous objects are simulated, but the partial differential equation modeling the object is too expensive to solve over the entire object, and a simpler, less expensive model would suffice in most of the object to reach the desired accuracy; air flow around an airplane is a typical example, where viscous effects are important close to the airplane, but can be neglected further away, see the early publication [START_REF] Dinh | On the coupling of viscous and inviscid models for incompressible fluid flows via domain decomposition[END_REF], and also [START_REF] Coclici | Analysis of a heterogeneous domain decomposition for compressible viscous flow[END_REF] and the references therein. An automatic approach for neglecting the diffusion in parts of the domain is the χ-formulation, see [START_REF] Lai | A defect equation approach for the coupling of subdomains in domain decomposition methods[END_REF] [START_REF] Brezzi | A self-adaptive formulation for the Euler/Navier-Stokes coupling[END_REF], and there are also techniques based on virtual control, originating in [START_REF]On the coupling of viscous and inviscid models for incompressible fluid flows via domain decomposition[END_REF], see [START_REF] Agoshkov | Optimal control in heterogeneous domain decomposition methods for advection-diffusion equations[END_REF] for the case with overlap, and [START_REF] Gervasio | Heterogeneous coupling by virtual control methods[END_REF] for the case without, and also [START_REF] Discacciati | heterogeneous mathematical models in fluid dynamics and associated solution algorithms[END_REF] for virtual control with variational coupling conditions. A third emerging area is the coupling of equations across dimensions, for example the blood flow in the artery can be modeled by a one dimensional model, but in the heart, it needs to be three dimensional, see for example [START_REF] Formaggia | On the coupling of 3d and 1d navier-stokes equations for flow problems in compliant vessels[END_REF]. All these techniques have become known in the domain decomposition community under the name heterogeneous domain decomposition methods, a terminology sparked by the review [START_REF] Quarteroni | Heterogeneous domain decomposition: principles, algorithms, applications[END_REF], and the literature has become vast in this field.

We are interested in this paper in the second situation, where the motivation for using different equations comes from the fact that we would like to use simpler, less expensive equations in areas of the domain where the full model is not needed, and we use as our guiding example the advection reaction diffusion equation. We are in principle interested in the fully viscous solution, but we would like to solve only an advection reaction equation for computational savings in part of the domain.

Coupling conditions for this type of problem have been developed in the seminal paper [START_REF] Gastaldi | On the coupling of hyperbolic and parabolic systems: analytical and numerical approach[END_REF], but with the first situation described above in mind, i.e. there is indeed a viscous and an inviscid physical domain, and the coupling conditions are obtained by a limiting process as the viscosity goes to zero, see also [START_REF] Gastaldi | On the coupling of two dimensional hyperbolic and elliptic equations : Analytical and numerical approach[END_REF], and [START_REF] Barbu | A class of singularly perturbed, coupled, evolutionary boundary value problems[END_REF][START_REF] Coclici | A heterogeneous domain decomposition for initial-boundary value problems with conservation laws and electromagnetic fields[END_REF] for an innovative correction layer, and [START_REF] Coclici | The coupling of hyperbolic and elliptic boundary value problems with variable coefficients[END_REF] for the steady case.

Dubach developed in his PhD thesis [START_REF] Dubach | Contribution à la Résolution des Équations fluides en domaine non borné[END_REF] coupling conditions based on absorbing boundary conditions, and such conditions have been used in order to define heterogeneous domain decomposition methods in [START_REF] Gander | Optimized Schwarz algorithms for coupling convection and convection-diffusion problems[END_REF]. A fundamental question however in the second situation described above is how far the solution obtained from the coupled problem is from the solution of the original, more expensive one on the entire domain. A first comparison of different transmission conditions focusing on this aspect appeared in [START_REF] Gander | Advection diffusion problems with pure advection approximation in subregions[END_REF]. In [START_REF] Gander | Viscous problems with inviscid approximations in subregions: a new approach based on operator factorization[END_REF], coupling conditions were developed for stationary advection reaction diffusion equations in one spatial dimensions, which lead to solutions of the coupled problem that can be exponentially close to the fully viscous solution, and rigorous error estimates are provided. The coupling conditions are based on the factorization of the differential operator, see also [START_REF] Lohéac | Parabolic approximations of the convectiondiffusion equation[END_REF], and the exact factorization can be used in this one dimensional steady case. We study in this paper time dependent advection reaction diffusion problems, where the exact factorization of the differential operator can not be used any more, due to the non-local nature of the factors, and new ideas are needed in order to obtain better coupling conditions than the classical ones developed for situation one, where the domains are really physically different.

We present in Section 2 our new factorization algorithm. In Section 3 and 4 we give a detailed error analysis of our algorithm, and prove asymptotic error estimates when the viscosity is becoming small. In Section 5 we present numerical experiments which show that our theoretical error estimates are sharp, and that the new factorization algorithm gives approximate solutions which are one order of magnitude more accurate in the viscous region than the best heterogeneous domain decomposition methods known from the literature.

2.

A new coupling algorithm based on factorization. We now explain how the factorization technique that led to coupling conditions of excellent quality for one dimensional problems in [START_REF] Gander | Viscous problems with inviscid approximations in subregions: a new approach based on operator factorization[END_REF] can be used to obtain a new coupling algorithm for evolution problems which we will call factorization algorithm.

Model problem. We consider the time dependent advection reaction diffusion equation

L ad u := ∂ t u -ν∂ 2 x u + a∂ x u + cu = f in (-L 1 , L 2 ) × (0, T ), B 1 u(-L 1 , •) = g 1 on (0, T ), B 2 u(L 2 , •) = g 2 on (0, T ), u(x, 0) = h in (-L 1 , L 2 ), (2.1)
where a is the velocity field, ν > 0 is the viscosity, and c > 0 is a reaction term. The B j , j = 1, 2 are suitable boundary operators, representing Dirichlet or absorbing boundary conditions. We consider the following choice of these operators, depending on the sign of the advection term a,

B 1 B 2 a > 0 Id ∂ t + a∂ x + c a < 0 Id Id (2.2)
In the case a > 0, the flow is given at the inflow boundary, and an absorbing boundary condition is prescribed at the outflow boundary. This can be compared to the situation of the tail of a wing, where the flow goes from the complicated model region into the simplified model region. For negative a, the flow is prescribed at the inflow and outflow boundary, which can be compared to the situation of the front of the wing, where the flow goes from the simplified model region into the complicated model region, and a boundary layer forms.

We assume that the initial data h is compactly supported in (-L 1 , 0). The forcing term f is compactly supported in (-L 1 , L 2 ) × (0, T ], and the boundary values g 1 and g 2 are compactly supported in (0, T ] . Regularity and compatibility conditions on the data need to be enforced to have a sufficiently regular solution, see Section 3.

2.2.

The new algorithm based on factorization. Using Nirenberg's factorization, we can factor the advection-diffusion operator into a product of two evolution operators in opposite x directions,

-ν∂ 2 x + a∂ x + c + ∂ t = -ν(∂ x + P + (∂ t )) (∂ x + P -(∂ t )) (mod C ∞ )
. This factorization has been used to design absorbing boundary conditions and paraxial equations for hyperbolic problems, see [START_REF] Bamberger | Higher order paraxial wave equation approximations in heterogeneous media[END_REF]. For parabolic problems, Nataf and coauthors [START_REF] Lohéac | Parabolic approximations of the convectiondiffusion equation[END_REF][START_REF] Nataf | Résolution de l'équation de convection-diffusion stationnaire par une factorisation parabolique[END_REF] computed approximations of u via a double sweep, and also obtained transmission conditions for Schwarz domain decomposition methods [START_REF] Nataf | Factorization of the convection-diffusion operator and the Schwarz algorithm[END_REF], which led to the new class of optimized Schwarz methods, see [START_REF] Gander | Optimized Schwarz methods[END_REF] for an overview. The same factorization can also be used to obtain incomplete LU preconditioners [START_REF] Gander | AILU: a preconditioner based on the analytic factorization of the elliptic operator[END_REF][START_REF] Gander | An incomplete lu preconditioner for problems in acoustics[END_REF], and is the underlying mathematical structure of the recently developed sweeping preconditioner [START_REF] Engquist | Sweeping preconditioner for the helmholtz equation: moving perfectly matched layers[END_REF]. We now use this factorization to define our new factorization algorithm: we define two subdomains, Ω 1 = (-L 1 , 0) and Ω 2 = (0, L 2 ), and want to couple the advection-diffusion equation in Ω 1 with an advection equation in Ω 2 , defined by the transport operator L a ≡ ∂ t + a∂ x + c. Our goal is to obtain a coupled solution which is as close as possible to the fully viscous solution of the original problem.

We start with the case a > 0, where according to (2.2) the exterior boundary condition is L a u(L 2 , •) = g 2 . Suppose there exists a decomposition L ad = L ma L a with L a a transport operator propagating to the right, and L ma a transport operator propagating to the left. The original problem

         L ma L a u = f in Ω × (0, T ), u(-L 1 , •) = g 1 on (0, T ), L a u(L 2 , •) = g 2 on (0, T ), u(•, 0) = h in Ω
can then be solved by introducing u ma := L a u, and solving the two problems

     L ma u ma = f in Ω 2 × (0, T ), u ma (L 2 , •) = g 2 on (0, T ), u ma (•, 0) = L a u(•, 0) in Ω 2 ,          L ma L a u ad = f in Ω 1 × (0, T ), u ad (-L 1 , •) = g 1 on (0, T ), L a u ad (0, •) = u ma (0, •) on (0, T ), u ad (•, 0) = h in Ω 1 , (2.3)
which leads to u ad = u |Ω1 . Unfortunately, the exact factorization L ad = L ma L a is expensive, but we can use an approximation with a remainder,

L ad = ν a 2 (L ma L a -R ) with R = (∂ t + c) 2 and L ma = ∂ t -a∂ x + c + a 2 ν . (2.4)
The viscous solution u satisfies L ma L a u = a 2 f /ν +Ru, and the algorithm corresponding to (2.3) is

         L ma u ma = a 2 ν f + R u in Ω 2 × (0, T ), u ma (L 2 , •) = g 2 on (0, T ), u ma (•, 0) = f (•, 0) + νd 2 x h in Ω 2 ,          L ad u ad = f in Ω 1 × (0, T ), u ad (-L 1 , •) = g 1 on (0, T ), L a u ad (0, •) = u ma (0, •) on (0, T ), u ad (•, 0) = h in Ω 1 .
Since u is unknown to evaluate the remainder, we approximate it by solving an advection equation, and our new factorization algorithm is

     L a u k a = f in Ω 2 × (0, T ), u k a (0, •) = u k-1 ad (0, •) on (0, T ), u k a (•, 0) = h in Ω 2 ,          L ma u k ma = a 2 ν f + R u k a in Ω 2 × (0, T ), u k ma (L 2 , •) = g 2 on (0, T ), u k ma (•, 0) = f (•, 0) + νd 2 x h in Ω 2 , (2.5)
           L ad u k ad = f in Ω 1 × (0, T ), u k ad (-L 1 , •) = g 1 on (0, T ), L a u k ad (0, •) = u k ma (0, •) on (0, T ), u k ad (•, 0) = h in Ω 1 ,
where we start with a given initial guess u 0 ad (0, •) = g 0 ad . We will prove well posedness of this algorithm in Section 3, and give precise error estimates when ν is small, which show that the new factorization algorithm gives one and a half orders of magnitude better solutions in the viscous subregion than the best other coupling algorithms from the literature.

When a < 0, we have the factorization with remainder in reverse order, L ad = ν a 2 (L a L ma -R ), and now the operator L a propagates to the left, and L ma to the right. The viscous solution u satisfies L a L ma u = a 2 f /ν + Ru, and introducing u a := L ma u, the algorithm corresponding to (2.3) is

       L a u a = a 2 ν f + R u in Ω 2 × (0, T ), u a (L 2 , •) = L ma u(L 2 , •) on (0, T ), u a (•, 0) = L ma u(•, 0) in Ω 2 ,          L ad u ad = f in Ω 1 × (0, T ), u ad (-L 1 , •) = g 1 on (0, T ), L ma u ad (0, •) = u a (0, •) on (0, T ), u ad (•, 0) = h in Ω 1 .
Since u is unknown to evaluate the remainder and the boundary conditions, we approximate it again by solving an advection equation, and our new factorization algorithm becomes

     L a u 1 a = f in Ω 2 × (0, T ), u 1 a (L 2 , •) = g 2 on (0, T ), u 1 a (•, 0) = h in Ω 2 ,          L a u 2 a = a 2 ν f + R u 1 a in Ω 2 × (0, T ), u 2 a (L 2 , •) = L ma u 1 a (L 2 , •) on (0, T ), u 2 a (•, 0) = L ma u 1 a (•, 0) in Ω 2 , (2.6)
         L ad u ad = f in Ω 1 × (0, T ), u ad (-L 1 , •) = g 1 on (0, T ), L ma u ad (0, •) = u 2 a (0, •) on (0, T ), u ad (•, 0) = h in Ω 1 ,
where one could also directly compute

L ma u 1 a (L 2 , •) = 2g ′ 2 + (2c + a 2 /ν)g 2 -f (L 2 , •) and L ma u 1 a (•, 0) = f (•, 0) -2ad x h + a 2 h/ν.
There is no iteration for a < 0 in the algorithm, because the boundary condition g 2 at x = L 2 in the first step can not be updated naturally from the viscous solution u ad in Ω 1 . We will study this algorithm in detail in Section 4, and show that it gives an order of magnitude better solutions in the viscous subregion than the other coupling algorithms from the literature.

Well-posedness results for advection reaction diffusion problems.

We work in the usual Sobolev spaces in time and space, H s (0, T ) and H s (Ω) for Ω ⊂ R, H s (Ω × (0, T )) in the hyperbolic case, and the anisotropic spaces H r,s (Ω × (0, T )) in the parabolic case. For clarity, we will add an index defining time or space in the Sobolev space, for instance H s t ≡ H s (0, T ). We introduce for any domain Ω ⊂ R the anisotropic Sobolev spaces (see [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]) H r,s (Ω × (0, T )) = L 2 (0, T ; H r (Ω)) ∩ H s (0, T ; L 2 (Ω)).

(2.7)

If u is in H r,s (Ω × (0, T )), then for any integer j and k, we have

∂ j ∂x j ∂ k ∂t k u ∈ H µ,ν (Ω × (0, T )), where µ r = ν s = 1 -( j r + k s ). (2.8)
We introduce the space V r,s of traces of functions in H r,s (Ω×(0, T )) for the half-space Ω = R -(and similarly for Ω = R + ). Denoting by f k the trace of the k-th derivative in time on the initial line, x ∈ R -, and by g j the trace of the j-th derivative in space on the boundary x = 0, t ∈ (0, T ), the trace space V r,s is defined by

V r,s := (f k , g j ) ∈ k<s-1 2 H p k (Ω) × j<r-1 2 H µj (0, T ), p k = r s (s -k -1 2 ), µ j = s r (r -j -1 2 ), ∂ k gj ∂t k (0) = ∂ j f k ∂x j (0), if j r + k s < 1 -1 2 ( 1 r + 1 s ), ∞ 0 | ∂ j f k ∂x j (σ s ) - ∂ k gj ∂t k (σ r )| 2 dσ σ < ∞, if j r + k s = 1 -1 2 ( 1 r + 1 s ) .
(2.9)

Theorem 2.1 ( [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]). For positive real numbers r and s such that 1-

1 2 ( 1 r + 1 s ) > 0, the trace map u → { ∂ k u ∂t k (x, 0)} k<s-1 2 , { ∂ j u ∂x j (0, t)} j<r-1 2 (2.10)
is defined and continuous from H r,s (Ω × (0, T )) onto V r,s . We start with well-posedness results for the advection equation, by stating a general result, applicable to L a in Ω or Ω 2 , and L ma in Ω 2 . To this end, we introduce O = (x 1 , x 2 ) and consider

L b v := ∂ t v + b ∂ x v + ηv = p in O × (0, T ). ( 2
= inf{t ≥ 0, s.t. x -at / ∈ Ō}. If b > 0, x -= x 1 and τ (x) = x-x1
b , and if b < 0, x -= x 2 and τ (x) = x-x2 b . Note that τ is a continuous function of x. We therefore equip (2.11) with initial and boundary conditions

v(•, 0) = h, v(x -, •) = g.
(2.12)

The following well-posedness result can be found in [START_REF] Métivier | Small viscosity and boundary layer methods[END_REF].

Theorem 2.2. If p ∈ L 2 (O × (0, T )), g ∈ L 2 (0, T ) and v 0 ∈ L 2 (O)
, then the transport problem (2.11,2.12) has a unique weak solution v ∈ L 2

x,t , given by (the characteristic function of ω in R 2 is denoted by

1 ω ) v(x, t) = h(x -bt)e -ηt 1 t<τ (x) + g(t -τ (x))e -ητ (x) 1 t>τ (x) + t (t-τ (x)) + p(x -b(t -s), s)e -η(t-s) ds.
(2.13)

If for some γ > 0 we have h ∈ H γ (O), g ∈ H γ (0, T ) and p ∈ H γ (O × (0, T ))
, with the compatibility conditions

d k t g(0) = ( k-1 j=0 (-M b ) j ∂ k-1-j t p)(x -, 0) + (-M b ) k h(x -) for 0 ≤ k ≤ γ -1, (2.14)
then v ∈ H γ (O×(0, T )) and v(x + , •) ∈ H γ (0, T ). Furthermore, we have for 0 ≤ k ≤ γ the estimates

η ∂ k t v 2 L 2 x,t + |b| ∂ k t v(x + , •) 2 L 2 t ≤ 1 η ∂ k t p 2 L 2 x,t + ∂ k t v(•, 0) 2 L 2 x + |b| d k t g 2 L 2 t . (2.15)
Similarly, we also use well-posedness results for the advection reaction diffusion equation

L ad u := ∂ t u -ν∂ 2 x u + a∂ x u + cu = f in O × (0, T ), B 1 u(x 1 , •) = g 1 on (0, T ), B 2 u(x 2 , •) = g 2 on (0, T ), u(x, 0) = h in O, (2.16)
with boundary operators according to (2.2). We define M ad to be the spatial part of the operator L ad , i.e.

L ad = ∂ t + M ad . Theorem 2.3. For γ > 0, let f ∈ H 2γ,γ (O × (0, T )), g 1 ∈ H γ+ 3 4 t , g 2 ∈ H γ+ 3 4 t
for negative a, and

g 2 ∈ H γ-1 4 t for positive a, h ∈ H 2γ+1 x (O), with the compatibility conditions for 0 ≤ k < γ -1 2 and 0 ≤ k ′ < γ -3 2 for negative a given by 1 ≤ j ≤ 2, d k t g j (0) = (-M ad ) k h(x j ) + ( k-1 j=0 (-M ad ) j ∂ k-1-j t f )(x j , 0), (2.17)
and for positive a the second compatibility condition is replaced by

d k ′ t g 2 (0) -ν k ′ -1 j=0 (-M ad ) j ∂ k ′ -1-j t ∂ 2 x f (x 2 , 0) -νd 2 x (-M ad ) k ′ h(x 2 ) = ∂ k ′ t f (x 2 , 0), (2.18)
then problem (2.16) has a unique solution u in H 2(γ+1),γ+1 (O × (0, T )).

Proof. Existence and regularity results are well-known for Dirichlet boundary conditions on both sides, see [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF][START_REF] Métivier | Small viscosity and boundary layer methods[END_REF], so we do not consider the case of negative advection further. In [START_REF] Métivier | Small viscosity and boundary layer methods[END_REF] more precise results with error bounds in ν for the hyperbolic equation (see Theorem 3.3) can be found. In the case where a > 0, due to the absorbing boundary, we need to modify the proof on the right boundary, and we use a Fourier transform in time. A weak solution is obtained by a variational formulation, like in [START_REF] Martin | An optimized Schwarz waveform relaxation method for unsteady convection diffusion equation[END_REF][START_REF] Bennequin | A homographic best approximation problem with application to optimized Schwarz waveform relaxation[END_REF] for instance. The regularity is obtained as follows: we first modify the boundary condition in (2.16) on the right at x = x 2 to Dirichlet data g2 ∈ H

γ+ 3 4 t
. Because of the compatibility conditions on the left, and imposing symmetric compatibility conditions on g2 on the right, there is a unique solution ũ ∈ H 2(γ+1),(γ+1) (O × (0, T )), see [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]. The difference v = u -ũ is solution of the homogeneous case of equation (3.1), but the boundary condition on the right becomes L a (u -ũ) = q 2 := g 2 -L a ũ.

By the regularity results above, q 2 is in H γ-1 4 t . To estimate v, we will make use of the Fourier transform. We extend all functions by 0 in R -, and smoothly into (T, +∞), and define

v(ω) = 1 2π R e -iω t v(t) dt.
Since the initial value vanishes, the equation is Fourier transformed in time to

L ad v := -ν∂ 2 x v + a∂ x v + (c + iω)v = 0 on O × C.
This is for each ω an ordinary differential equation, with characteristic roots

λ + (ω) = 1 2ν (a + a 2 + 4ν(c + iω)), λ -(ω) = 1 2ν (a -a 2 + 4ν(c + iω)), (2.19)
with Re(λ + ) > 0 and Re(λ -) < 0. The general solution is

v(x, ω) = ℓ + (ω)e λ+x + ℓ -(ω)e λ-x .
Using the boundary conditions, we then get the solution

v(x, ω) = q2 (ω) e λ+(x-x1) -e λ-(x-x1) νλ 2 + e λ+(x2-x1) -νλ 2 -e λ-(x2-x1) , (2.20)
where we have used the relation c + iω + aλ ± = ν(λ ± ) 2 . The value at x = x 2 can be equivalently written as

v(x 2 , ω) = q2 (ω) e -(λ+-λ-)(x2-x1) -1 νλ 2 + ( λ- λ+ ) 2 e -(λ+-λ-)(x2-x1) -1 . (2.21)
In order to estimate the regularity of v(x 2 , •), we need to estimate the multiplicative factor on the right for large ω. We can see from (2.19)

that λ + (ω) ∼ -λ -(ω) ∼ iω ν . Therefore |v(x 2 , ω)| ∼ q2(ω) νλ 2 + ∼ q2(ω) ω . Since q 2 ∈ H γ-1 4 t , we conclude that v(x 2 , •) ∈ H γ+ 3 4 t
. Then v is solution of the advection-diffusion equation with Dirichlet boundary conditions, and the data has sufficient regularity to conclude.

3. Properties of the factorization algorithm for positive advection. We consider the advection-diffusion equation in Ω = (-L 1 , L 2 ) with Dirichlet boundary condition on the left, and absorbing boundary condition given by the transport operator on the right (see [START_REF] Halpern | Artificial boundary conditions for the linear advection diffusion equation[END_REF]),

L ad u := ∂ t u -ν∂ 2 x u + a∂ x u + cu = f in Ω × (0, T ), u(-L 1 , •) = g 1 on (0, T ), L a u(L 2 , •) = g 2 on (0, T ), u(•, 0) = h in Ω. (3.1)
We suppose in this section that f and g 1 are compactly supported in (0, T ], that h is compactly supported in Ω 1 = (-L 1 , 0), and that for each t the function f (•, t) is compactly supported in Ω. We further assume that the boundary condition at L 2 is absorbing, that is g 2 = 0. Therefore the compatibility conditions are satisfied to any order on both ends of the interval Ω, and for f ∈ H

9 2 , 9 4 (Ω × (0, T )), h ∈ H 11 2
x , and 13 4 (Ω × (0, T )). 3.1. Well-posedness. The remainder R for computing u k ma in the new factorization algorithm (2.5) contains two time derivatives, which lead to an important loss of regularity at each iteration. We will however see that the error order in ν can not be improved further after two iterations, and hence we only study the first two iterations in detail. We start with the well-posedness of the algorithm. Algorithm (2.5) starts with an initial guess g 0 ad as boundary condition for u a . We assume that g 0 ad ∈ H 4 t and is compactly supported in (0, T ]. Using that f ∈ H 4 (Ω 2 × (0, T )), that h vanishes in Ω 2 , and that the compatibility conditions at x = 0, t = 0 are satisfied, the solution of

g 1 ∈ H 3 t , u is defined in H 13 2 ,
L a u 1 a = f in Ω 2 , u 1 a (0, •) = g 0 ad , u 1 a (•, 0) = 0 satisfies u 1 a ∈ H 4 (Ω 2 × (0, T )
). The right hand side for the modified advection equation in (2.5) is then

f 1 ma = a 2 ν f + R u 1 a ∈ H 2 (Ω 2 × (0, T ))
, and solving

L ma u 1 ma = f 1 ma in Ω 2 , u 1 ma (L 2 , •) = 0, u 1 ma (•, 0) =
0, the compatibility conditions at x = L 2 are again satisfied to any order, which implies that u 1 ma ∈ H 2 (Ω 2 × (0, T )) and u 1 ma (0, •) ∈ H 2 (0, T ). The latter then becomes the right boundary data for the advection diffusion problem in Ω 1 ,

L ad u 1 ad = f in Ω 1 , u 1 ad (-L 1 , •) = g 1 , L a u 1 ad (0, •) = u 1 ma (0, •), u 1 ad (•, 0) = h.
We have seen already that the compatibility conditions on the left are satisfied, and on the right, at the corner (0, 0), with the regularity of u 1 ma , the condition u 1 ma (0, 0)νd 2

x h(0) = f (0, 0) holds, since both sides of this equality vanish. Since f ∈ H

9 2 , 9 4 (Ω × (0, T )), h ∈ H 11 2
x , g 1 ∈ H 3 t , and u 1 ma (0, •) ∈ H 2 t , we obtain u 1 ad ∈ H 13 2 , 13 4 (Ω 1 × (0, T )) and u 1 ad (0, •) ∈ H 3 t , and at the corner (0, 0), we have for g 1 ad := u 1 ad (0, •)

g 1 ad (0) = h(0), d t g 1 ad (0) + M ad h(0) = f (0, 0), d 2 t g 1 ad (0) -M 2 ad h(0) = ∂ t f (0, 0) -M ad f (0, 0). (3.2)
We now start the second iteration with the computation of u 2 a , using u 1 a (0,

•) = g 1 ad = u 1 ad (0, •) ∈ H 3 t .
Since h is compactly supported in Ω 1 , M p ad h(0) = M p a h(0) = 0, and (3.2) are appropriate compatibility conditions to compute u 2 a ∈ H 3 (Ω 2 × (0, T )). We define the new right hand side

f 2 ma = a 2 ν f + R u 2 a ∈ H 1 (Ω 2 × (0, T )),
and compute the solution of

L ma u 2 ma = f 2 ma in Ω 2 , u 2 ma (L 2 , •) = 0, u 2 ma (•, 0) = 0.
We thus obtain u 2 ma ∈ H 1 (Ω 2 × (0, T )) and u 2 ma (0, •) ∈ H 1 (0, T ). The last step of the second iteration is to compute u 2 ad solution of

L ad u 2 ad = f in Ω 1 , u 1 ad (-L 1 , •) = g 1 , L a u 1 ad (0, •) = u 2 ma (0, •), u 1 ad (•, 0) = h,
and we need only to satisfy a compatibility condition on the left, which implies that 4 (Ω × (0, T )).

u 2 ad ∈ H 9 2 , 9

3.2.

Error estimates for the factorization algorithm. We present now asymptotic error estimates for small viscosity ν when u, the viscous solution of (3.1), is approximated by (u k ad , u k a ), the solution obtained by our new factorization algorithm (2.5). We define the error quantities e k a := u k a -u, e k ad := u k ad -u, e k ma := u k ma -L a u := u k ma -u ma , and suppose that all our data is C ∞ in all variables. The error equations are

   L a e k a = -ν∂ 2 x u in Ω 2 , e k a (0, •) = e k-1 ad (0, •), e k a (•, 0) = 0,    L ma e k ma = R e k a in Ω 2 , e k ma (L 2 , •) = 0, e k ma (•, 0) = 0,        L ad e k ad = 0 in Ω 1 , e k ad (-L 1 , •) = 0, L a e k ad (0, •) = e k ma (0, •), e k ad (•, 0) = 0, (3.3
) with e 0 ad (0, •) := g 0 ad -u(0, •). We need more precise estimates than those provided by Theorems 2.2 and 2.3. First, we state precisely the initial conditions for all the equations involved: the parabolic problems in Ω and Ω 1 will use

∂ k t u(•, 0) = (-M ad ) k h and ∂ k t u p ad (•, 0) = (-M ad ) k h, the forward hyperbolic problem in Ω 2 uses ∂ k t u p a (•, 0) = (-M a ) k h = 0,
and the backward hyperbolic problem in Ω 2 uses

∂ k t u p ma (•, 0) = k-1 j=0 (-M ma ) j ∂ k-1-j t f p ma (•, 0) + (-M ma ) k u ma,0 = 0.
The solution of the exact backward hyperbolic problem in Ω 2 has as initial condition

∂ k t u ma (•, 0) = k-1 j=0 (-M ma ) j ∂ k-1-j t f ma (•, 0) + (-M ma ) k u ma,0 = 0,
from which we infer the initial values for the errors such that

∂ k t e p ad (•, 0) = 0 in Ω 1 , ∂ k t e p a (•, 0) = 0 and ∂ k t e p ma (•, 0) = 0 in Ω 2 . (3.4)
We start with estimates for the solution of the advection-diffusion equation (3.1) with vanishing initial data and vanishing boundary data g 1 . A first lemma gives results for the problem with g 2 = 0, based on energy estimates, and a second lemma gives estimates where only the right-hand side f is non-zero.

Lemma 3.1. Suppose that a > 0, and that h vanishes identically in Ω, g 1 and g 2 vanish on (0, T ), and that f is in C ∞ 0 (Ω × (0, T ]). Then there is a positive constant C such that for any ν > 0, and any k ≤ γ, the solution u 1 of (3.1) satisfies the estimates

∂ k t u 1 2 L 2 x,t + ∂ k t u 1 (L 2 , •) 2 L 2 t + ν ∂ k t ∂ x u 1 2 L 2 x,t ≤ C ∂ k t f 2 L 2 x,t , (3.5)
ν ∂ k t ∂ 2 x u 1 L 2 x,t ≤ ∂ k t f L 2 x,t . (3.6)
Proof. Since the compatibility conditions are satisfied, u 1 is in H ∞ , and the initial value of ∂ k t u 1 vanishes as well. We start with k = 0: multiplying the equation by u 1 and integrating over Ω, taking into account that u 1 vanishes at -L 1 gives 1 2

d dt u 1 (•, t) 2 L 2 x + c u 1 (•, t) 2 L 2 x + ν ∂ x u 1 (•, t) 2 L 2 x + a 2 u 2 1 (L 2 , t) -ν(u 1 ∂ x u 1 )(L 2 , t) = Ω f (x, t)u 1 (x, t)dx.
Using the boundary condition at L 2 yields 1 2

d dt ( u 1 (•, t) 2 L 2 x + ν a u 2 1 (L 2 , t)) + c u 1 (•, t) 2 L 2 x + ν ∂ x u 1 (•, t) 2 L 2 x + ( a 2 + νc a )u 2 1 (L 2 , t) = Ω f (x, t)u 1 (x, t)dx,
and by Cauchy-Schwarz and Young's inequality we obtain

Ω f (x, t)u 1 (x, t)dx ≤ u 1 (•, t) L 2 x f (•, t) L 2 x ≤ c 2 u 1 (•, t) 2 L 2 x + 1 2c f (•, t) 2 L 2 x .
Integrating in time over (0, T ), and dropping the first term which is positive, we obtain, since the initial data vanishes, the inequality

c u 1 2 L 2 x,t + 2ν ∂ x u 1 2 L 2 x,t + a u 1 (L 2 , •) 2 L 2 t ≤ 1 c f 2 L 2 x,t
, which proves (3.5) for k = 0. To prove (3.6), we multiply the equation by

-∂ 2 x u 1 and integrate in x, ν ∂ 2 x u 1 (•, t) 2 L 2 x -(∂ t u 1 (•, t), ∂ 2 x u 1 (•, t)) -a(∂ x u 1 (•, t), ∂ 2 x u 1 (•, t)) -c(u 1 (•, t), ∂ 2 x u 1 (•, t)) = -(f (•, t), ∂ 2 x u 1 (•, t)).
An integration by parts leads to

ν ∂ 2 x u 1 (•, t) 2 L 2 x + 1 2 d dt ∂ x u 1 (•, t) 2 L 2 x + c ∂ x u 1 (•, t) 2 L 2 x -[∂ x u 1 (•, t)(∂ t u 1 + a 2 ∂ x u 1 + cu 1 )(•, t)] L2 -L1 = -(f (•, t), ∂ 2 x u 1 (•, t)).
By the boundary conditions, the boundary terms become

a(∂ x u 1 ) 2 (-L 1 , t) + 1 a (∂ t u 1 + cu 1 ) 2 (L 2 , t) > 0.
We can now integrate in time and use Cauchy-Schwarz and Young's inequality to obtain

ν ∂ 2 x u 1 2 L 2 x,t ≤ 1 ν f 2 L 2 x,t .
The estimates with the time derivative are obtained by applying the equation to ∂ k t u.

Lemma 3.2. Assume that a > 0. Then there are constants ν > 0 and C > 0 such that for ν ≤ ν, and for any g 2 ∈ C ∞ 0 ((0, T ]), the solution u 2 of (3.1) with zero data h, f and g 1 satisfies for all k ≤ γ the inequalities

∀x ∈ [-L 1 , L 2 ], ∂ k t u 2 (x, •) L 2 t ≤ Cν g 2 H k t , (3.7)
∂ k t u 2 L 2 x,t ≤ Cν 3 2 g 2 H k t , ∂ x ∂ k t u 2 L 2 x,t ≤ Cν 1 2 g 2 H k t , ∂ 2 x ∂ k t u 2 L 2 x,t ≤ Cν -1 2 g 2 H k t .
Proof. We use a Fourier transform argument as in the proof of Theorem 2.3, and rewrite (2.20) as û2 (x, ω) = ĝ2 (ω) e λ+(x-L2)

e -(λ+-λ-)(x+L1) -1

νλ 2 + ( λ- λ+ ) 2 e -(λ+-λ-)(L2+L1) -1 . (3.8)
Now we search for estimates in ν that are uniform in ω. We have for the roots λ ± the estimates

|λ -/λ + | < 1, Re(λ + -λ -) ≥ a/ν, |λ + | ≥ a/ν.
The numerator in (3.8) is bounded by 2. A lower bound for the denominator is obtained by writing |νλ 2 + | ≥ a 2 ν together with

1 - λ - λ + 2 e -(λ+-λ-)(L2+L1) ≥ 1 - λ - λ + 2 e -Re(λ+-λ-)(L2+L1) ≥ 1 -e -a ν (L2+L1) .
Inserting these estimates into (3.8) gives

|û 2 (x, ω)| ≤ 2ν a 2 |ĝ 2 (ω)|e Reλ+(x-L2) 1 1 -e -a ν (L2+L1) .
Since 1/|1 -µ| < 2 for µ < 1/2, for ν sufficiently small so that e -a ν (L2+L1) < 1/2, we have for any ω,

|û 2 (x, ω)| ≤ 4ν a 2 |ĝ 2 (ω)|e Reλ+(x-L2) ≤ Cν |ĝ 2 (ω)|. (3.9)
By Parseval's identity, we obtain

u 2 (x, •) L 2 (R+) ≤ Cν g 2 L 2 (R+) . (3.10)
Modifying now g 2 to vanish in [T + ǫ, ∞), the solution in (0, T ) remains unaffected by causality and for any positive ǫ,

u 2 (x, •) L 2 (0,T ) ≤ Cν g 2 L 2 (0,T +ǫ) .
Since ǫ is arbitrary, we conclude that

u 2 (x, •) L 2 (0,T ) ≤ Cν g 2 L 2 (0,T ) .
From (3.9) we also obtain for all ω

û2 (•, ω) 2 L 2 x ≤ Cν 2 Reλ + |ĝ 2 (ω)| 2 ≤ Cν 3 |ĝ 2 (ω)| 2 .
We thus obtain

u 2 L 2 x,t ≤ Cν 3 2 g 2 L 2 t .
For the derivative in space, we compute

∂ x û2 (x, ω) = ĝ2 (ω) λ+e λ + (x+L 1 ) -λ-e λ -(x+L 1 ) νλ 2 + e λ + (L 2 +L 1 ) -νλ 2 -e λ -(L 2 +L 1 ) = ĝ2 (ω) e λ+(x-L2) λ - λ + e -(λ + -λ -)(x+L 1 ) -1 νλ+ ( λ - λ +
) 2 e -(λ + -λ -)(L 2 +L 1 ) -1

.

For small ν, we therefore get as before

|∂ x û2 (x, ω)| ≤ 4 ν|λ + | |ĝ 2 (ω)| e Reλ+(x-L2) ≤ 4 a |ĝ 2 (ω)| e Reλ+(x-L2) , which gives ∂ x û2 2 L 2 (Ω×R) ≤ Cν ĝ2 2 L 2 (R)
, and using the same arguments as before,

∂ x u 2 L 2 x,t ≤ Cν 1 2 g 2 L 2 t .
In the same way we find that

∂ 2 x u 2 L 2 x,t ≤ Cν -1 2 g 2 L 2 t .
Theorem 3.3. Let a > 0. Then there exist positive constants C and ν such that for any ν ≤ ν, and for any set of data h ∈ C ∞ 0 (Ω), f ∈ C ∞ 0 (Ω×(0, T ]), g 1 ∈ C ∞ 0 ((0, T ]) and g 2 ≡ 0, if U is the solution of the transport equation

L a U = f in Ω × (0, T ), U (•, 0) = h U (-L 1 , •) = g 1 , (3.11)
then the solution u of the advection-diffusion equation (3.1) satisfies the estimate

∂ k t (u -U ) 2 L 2 x,t + ∂ k t (u -U )(L 2 , •) 2 L 2 t + ν ∂ k t ∂ x (u -U ) 2 L 2 x,t (3.12) + ν 2 ∂ k t ∂ 2 x (u -U ) 2 L 2 x,t ≤ Cν 2 ∂ k t ∂ 2 x U 2 L 2 x,t
.

Hence u also satisfies the estimate

∂ k t u 2 L 2 x,t + ∂ k t u(L 2 , •) 2 L 2 t + ∂ k t ∂ x u 2 L 2 x,t (3.13) + ∂ k t ∂ 2 x u 2 L 2 x,t ≤ C( f 2 H k+2 x,t + h 2 H k+2 x + g 1 2 H k+2 t ).
Proof. Since the data is compactly supported, the compatibility conditions are automatically satisfied, so that the solutions of the parabolic and hyperbolic equations are in C ∞ 0 (Ω × (0, T ]). The estimates (3.12) follow directly from Lemma 3.1, using that u -U is solution of the advection-diffusion equation in Ω with right hand side ν∂ 2

x U , and zero initial and boundary conditions on the left. The boundary condition on the right also vanishes,

L a (u -U ) = -f (L 2 , •) = 0, since f is compactly supported in Ω. We define now B k = f 2 H k x,t + h 2 H k x + g 1 2 H k t
, and use for U the hyperbolic estimates (2.15) in O = Ω,

∂ k t U 2 L 2 x,t + ∂ k t U (L 2 , •) 2 L 2 t ≤ CB k .
Next, from the advection equation, we deduce that

a∂ x U = -(∂ t + c)U + f, and a 2 ∂ 2 x U = (∂ t + c) 2 U + (a∂ x -(∂ t + c))f, so that ∂ k t ∂ x U 2 L 2 x,t ≤ CB k+1 , ∂ k t ∂ 2 x U 2 L 2 x,t ≤ CB k+2 ,
and from (3.12) we obtain

∂ k t (u -U ) 2 L 2 x,t + ∂ k t (u -U )(L 2 , •) 2 L 2 t ≤ Cν 2 B k+2 , ∂ k t ∂ x (u -U ) 2 L 2 x,t ≤ CνB k+2 , ∂ k t ∂ 2 x (u -U ) 2 L 2 x,t ≤ CB k+2 . Writing u = u -U + U gives ∂ k t u 2 L 2 x,t + ∂ k t u(L 2 , •) 2 L 2 t ≤ C(B k + ν 2 B k+2 ), ∂ k t ∂ x u 2 L 2 x,t ≤ C(B k+1 + νB k+2 ), ∂ k t ∂ 2 x u 2 L 2 x,t ≤ CB k+2 .
Therefore there is a new constant C and ν such that for ν ≤ ν, (3.13) holds. We now present an improved estimate for the solution of the modified advection problem in Ω 2 .

Theorem 3.4. Let a > 0. Then there exist positive constants C and ν such that, for ν ≤ ν, and for any set p compactly supported in Ω 2 × (0, T ], the solution v of the initial boundary value problem with modified advection   

L ma v = p in (0, L 2 ) × (0, T ), v(L 2 , •) = 0 on (0, T ), v(•, 0) = 0 in (0, L 2 )
satisfies the estimate

∂ k t v(0, •) 2 L 2 t ≤ Cν 2 ∂ k t p(0, •) 2 L 2 t + e -2 aL 2 ν ∂ k t p(L 2 , •) 2 L 2 t + ν ∂ k t p 2 H 1 x,t . (3.14)
Proof. We first extend p by 0 on (T, +∞). As in Theorem 2.2, v can be obtained using the method of characteristics,

v(0, t) = t max(t- L 2 a ,0)
p(a(t -σ), σ)e -c(t-σ) dσ, where c = c + a 2 /ν. Integrating by parts, and denoting by d t the characteristic derivative, i.e. d t φ = ∂ t φ -a∂ x φ, we obtain v(0, t) = 1 c p(0, t)

I - 0 for t < L 2 /a, p(L 2 , t -L2 a )e -c L 2 a for t > L 2 /a II - t max(t- L 2 a ,0) d t p(a(t -σ), σ)e -c(t-σ) dσ III .
The norm of the first term is

I L 2 t = p(0, •) L 2 t
, and the norm of the second term can be estimated as

II 2 L 2 t = +∞ L 2 a p 2 (L 2 , t - L 2 a )e -2c L 2 a dt ≤ e -2cL2/a p(L 2 , •) 2 L 2 t .
For the norm of the third term, we get

III 2 L 2 t = L 2 a 0 t 0 d t p(a(t -σ), σ)e -c(t-σ) dσ 2 dt + +∞ L 2 a t t- L 2 a d t p(a(t -σ), σ)e -c(t-σ) dσ 2 dt,
and using the Cauchy-Schwarz inequality, we obtain

III 2 L 2 t ≤ 1 2c L 2 a 0 t 0 (d t p) 2 (a(t -σ), σ) dσ dt + +∞ L 2 a t t- L 2 a (d t p) 2 (a(t -σ), σ) dσ dt ≤ 1 2c +∞ 0 L2 0 (d t p) 2 (x, t)dx dt ≤ 1 2c p 2 H 1 x,t , (3.15)
which finally leads to the estimate

v(0, •) 2 L 2 t ≤ C c2 p(0, •) 2 L 2 t + e -2c L 2 a p(L 2 , •) 2 L 2 t + 1 2c p 2 H 1 x,t . Since c ≥ a 2 ν , we get v(0, •) 2 L 2 t ≤ Cν 2 p(0, •) 2 L 2 t + e -2 aL 2 ν p(L 2 , •) 2 L 2 t + ν p 2 H 1 x,t
on the enlarged time interval (0, +∞). Using that the extension is vanishing for t ≥ T + ǫ gives the estimate (3.14) for k = 0 for any ǫ, and thus on (0, T ). Applying (3.14) to ∂ k t v gives then the general result. Theorem 3.5. Assume that a > 0, and let

B k := f 2 H k x,t + h 2 H k x + g 1 2 H k t .
Then there exist positive constants C and ν such that, for any data

h ∈ C ∞ 0 (Ω 1 ), f ∈ C ∞ 0 (Ω × (0, T ]), g 1 ∈ C ∞ 0 ((0, T ]
) and g 2 ≡ 0, and for any initial guess g 0 ad ∈ C ∞ 0 ((0, T ]) and any ν ≤ ν, the approximation from the new algorithm (2.5) satisfies the error bounds

u -u 1 a 2 L 2 x,t ≤ C(B 2 + g 0 ad 2 L 2 t ), u -u 1 ad 2 L 2 x,t ≤ Cν 5 (B 5 + g 0 ad 2 H 3 t ) (3.16) u -u 2 a 2 L 2 x,t ≤ Cν 2 (B 5 + ν 2 g 0 ad 2 H 3 t ), u -u 2 ad 2 L 2 x,t ≤ Cν 8 (B 8 + ν g 0 ad 2 H 6 t
).

(3.17)

Proof. We start with the proof for the first iteration of the new algorithm (2.5), with the initial guess g 0 ad , which gives e 0 ad = g 0 ad -u(0, •) in the algorithm (3.3) satisfied by the errors:

Advection:

The error e 1 a is solution of an advection equation in Ω 2 with right hand side -ν∂ 2

x u. Since the initial conditions vanish as described in (3.4), the hyperbolic estimate in Theorem 2.2 gives

∂ k t e 1 a 2 L 2 x,t + ∂ k t e 1 a (L 2 , •) 2 L 2 t ≤ C(ν 2 ∂ k t ∂ 2 x u 2 L 2 x,t + e 0 ad 2 H k t ).
We bound e 0 ad H k t by g 0

ad H k t + u(0, •) H k t . From (3.13), for small ν, ∂ k t u(0, •) 2 L 2
x,t

(by the trace theorem) and

∂ k t ∂ 2 x u 2 L 2
x,t are bounded by CB k+2 , which gives

∂ k t e 1 a 2 L 2 x,t + ∂ k t e 1 a (L 2 , •) 2 L 2 t ≤ C(ν 2 B k+2 + e 0 ad 2 H k t ) ≤ C(ν 2 B k+2 + g 0 ad 2 H k t +B k+2 ).
(3.18) This equation gives for k 0 and small ν the first estimate in (3.16).

Modified advection: Let p 1 a (0, •) := (∂ t +c) 2 e 1 a (0, •) = (∂ t +c) 2 e 0 ad . We estimate k t e 1 ma at x = 0 using Theorem 3.4,

∂ k t e 1 ma (0, •) 2 L 2 t ≤ Cν 2 ( ∂ k t p 1 a (0, •) 2 L 2 t + e -2 aL 2 ν ∂ k t p 1 a (L 2 , •) 2 L 2 t + ν ∂ k t p 1 a 2 H 1 x,t ) ≤ Cν 2 ( e 0 ad 2 H k+2 t + e -2 aL 2 ν ∂ k t p 1 a (L 2 , •) 2 L 2 t + ν ∂ k t p 1 a 2 H 1 x,t
). (3.19) For the last term on the right, we obtain

∂ k t p 1 a 2 H 1 x,t = ∂ k t p 1 a 2 L 2 x,t + ∂ x ∂ k t p 1 a 2 L 2 x,t + ∂ k+1 t p 1 a 2 L 2 x,t ≤ C( p 1 a 2 H k+1 (0,T ;L 2 x ) + ∂ x p 1 a 2 H k (0,T ;L 2 x ) ) ≤ C( e 1 a 2 H k+3 (0,T ;L 2 x ) + ∂ x e 1 a 2 H k+2 (0,T ;L 2 x ) ). Since L a e 1 a = -ν∂ 2 x u, we get ∂ x e 1 a = -1 a (ν∂ 2
x u+(∂ t +c)e 1 a ), and hence ∂ x e 1 a H k+2 (0,T ;L 2

x ) ≤ C(ν ∂ 2

x u H k+2 (0,T ;L 2 x ) + e 1 a H k+3 (0,T ;L 2

x ) ). Using (3.18), we finally obtain

∂ k t p 1 a 2 H 1 x,t ≤ C(ν 2 B k+5 + e 0 ad 2 H k+3 t + ν 2 ∂ 2 x u 2 H k+2 (0,T ;L 2 x )
). We use now (3.13) which gives ∂ 2

x u 2 H k+2 (0,T ;L 2

x ) ≤ CB k+4 and return to (3.19). From the hyperbolic estimate (3.18), we see that

∂ k t e 1 a (L 2 , •) 2 L 2 t and ∂ k t e 1 a 2 L 2
x,t are bounded by the same quantity, and hence we can also use the same bound for e -2 aL 2 ν

∂ k t p 1 a (L 2 , •) 2 L 2 t and ν ∂ k t p 1 a 2 H 1 x,t , ∂ k t e 1 ma (0, •) 2 L 2 t ≤C(ν 2 e 0 ad 2 H k+3 t + ν 5 B k+5 ). (3.20)
Advection-diffusion: e 1 ad is solution of the advection-diffusion equation with non-zero data only on the right. Therefore, applying Lemma 3.2 in Ω 1 with g 2 = e 1 ma (0, •), we obtain We now prove the error estimates for the second iteration: Advection: We again use the hyperbolic estimates for e 2 a . Since the initial values are also vanishing, we obtain as in (3.18) the estimate

e 1 ad 2 L 2 x,t ≤ Cν 3 e 1 ma (0, •)
∂ k t e 2 a 2 L 2 x,t + ∂ k t e 2 a (L 2 , •) 2 L 2 t ≤ C(ν 2 B k+2 + e 1 ad (0, •) 2 H k t ). (3.24)
Inserting (3.22) we get

∂ k t e 2 a (L 2 , •) 2 L 2 t + ∂ k t e 2 a 2 L 2 x,t ≤ C(ν 2 B k+2 + ν 4 (ν 3 B k+5 + e 0 ad 2 H k+3 t )) ≤ Cν 2 (B k+5 + ν 2 e 0 ad 2 H k+3 t ) (3.25) ≤ Cν 2 (B k+5 + ν 2 g 0 ad 2 H k+3 t + ν 2 B k+5 ).
The last estimate with k = 0 gives the first result in (3.17).

Modified advection: Defining p 2 a := R e 2 a = (∂ t + c) 2 e 2 a , we obtain using (3.25)

∂ k t p 2 a (L 2 , •) 2 L 2 t + ∂ k t p 2 a 2 L 2 x,t ≤ Cν 2 (B k+7 + e 0 ad 2 H k+5 t ). (3.26)
We estimate e 2 ma at x = 0 by Theorem 3.4,

e 2 ma (0, •) 2 L 2 t ≤ Cν 2 ( p 2 a (0, •) 2 L 2 t + e -2 aL 2 ν p 2 a (L 2 , •) 2 L 2 t + ν p 2 a 2 H 1 x,t
). (3.27) As in the first step, the term at the boundary x = L 2 is absorbed in the volume term, and p 2 a (0, •) = R e 1 ad (0, •), which can be estimated by (3.23). To estimate the term p 2 a H 1 x,t , we proceed as in the first iteration, to obtain ).

p 2 a 2 H 1 x,t ≤ C( e 2 a 2 H 3 (0,T ;L 2 x ) + ν 2 ∂ 2 x u 2 H 2 (0,T ;L 2 x ) ) ≤ C( e 2 a 2 H 3 (0,T ;L 2 x ) + ν 2 B 4 ) ≤ Cν 2 (B 8 + ν 2 e 0 ad 2 H 6 t ). ( 3
4. Properties of the factorization algorithm for negative advection. We consider now the advection-diffusion equation for a < 0 in Ω = (-L 1 , L 2 ) with Dirichlet boundary conditions on both sides,

L ad u := ∂ t u -ν∂ xx u + a∂ x u + cu = f in Ω × (0, T ), u(-L 1 , •) = g 1 in (0, T ), u(L 2 , •) = g 2 in (0, T ), u(•, 0) = h in Ω. (4.1)
We suppose again that f and (g 1 , g 2 ) are compactly supported in (0, T ], and that h is compactly supported in Ω 1 = (-L 1 , 0), and that for each t the function f (•, t) is compactly supported in Ω.

4.1. Well-posedness of the factorization algorithm. For u 1 a , suppose that f ∈ H 3+ 3 4 (Ω × (0, T )), g 2 ∈ H 3+ 3 4 (0, T ), h ∈ H 3+ 3 4 (Ω), and that the compatibility conditions (2.14) are satisfied. Then we have a unique solution u 1 a in W 3+ 3 4 (Ω 2 × (0, T )). For u 2 a , using the previous result, we have f -R u 1 a ∈ H 1+ 3 4 (Ω 2 × (0, T )), and 4 (Ω 2 × (0, T )). Therefore the traces at x = L 2 and t = 0 are in H 1+ 4 (Ω 1 × (0, T )). Finally for u, using the regularity assumptions above, u ∈ H 2(γ+1),γ+1 (Ω × (0, T )) with γ = 11 8 . 4.2. Error estimates for the factorization algorithm. We need a further lemma in order to obtain our asymptotic estimates.

L ma u 1 a ∈ H 2+ 3
Lemma 4.1. Suppose a < 0, and let g ∈ L 2 (0, T ). Then there exists a constant C > 0, such that for all ν > 0 the solution v of

L ad v = 0 in Ω 1 × (0, T ), v(-L 1 , •) = 0 on (0, T ), (∂ t -a∂ x + ( a 2 ν + c))v(0, •) = g on (0, T ), v(•, 0) = 0 in Ω 1 , (4.2)
satisfies the a priori estimate

v 2 L 2 x,t ≤ Cν 2 g 2 L 2 t .
Proof. Multiplying the equation by v, integrating on (-L 1 , 0) and using the boundary condition at

x = -L 1 yields 1 2 d dt v(•, t) 2 L 2 x - |a| 2 v 2 (0, t) + ν ∂ x v(•, t) 2 L 2 x -ν∂ x v(0, t)v(0, t) + c v(•, t) 2 L 2 x = 0. Inserting the boundary condition at x = 0 we obtain 1 2 d dt ( v(•, t) 2 L 2 x + ν |a| v 2 (0, t)) + ( |a| 2 + νc |a| )v 2 (0, t) +ν ∂ x v(•, t) 2 L 2 x + c v(•, t) 2 L 2 x = ν |a| g(t)v(0, t).
Using the inequality ν|g(t)v(0, t)|/|a| ≤ ν 2 |a| 3 g 2 (t) + |a| 4 v 2 (0, t) and integrating on the time interval (0, T ) gives for all t ∈ (0, T )

c v 2 L 2 x,t ≤ ν 2 |a| 3 t 0 g 2 (τ ) dτ.
We can now prove our main theorem for negative advection. Theorem 4.2. Suppose a < 0. Then there are positive constants C and ν such that for any h ∈ C ∞ 0 (Ω 1 ), f ∈ C ∞ 0 (Ω × (0, T ]), g 1 , g 2 ∈ C ∞ 0 ((0, T ]), and for any ν ≤ ν, the solution obtained by the new factorization algorithm (2.6) satisfies the estimates

u -u 1 a L 2 x,t ≤ Cν ∂ 2 x u L 2 x,t , (4.3
) u -u ad L 2 x,t ≤ Cν 2 ( u H 2,2 x,t + f (•, 0) H 2 x + h H 4 x + ∂ 2 x u(L 2 , •) L 2 t ), (4.4
)

which implies that u -u 1 a L 2 x,t ν, u -u ad L 2 x,t ν 2 .
Proof. We define the errors e 1 a := u 1 a -u, e ad := u ad -u, and e 2 a := u 2 a -L ma u.

Since L ma u(•, 0) = f (•, 0) -2ad x h + a 2 h/ν + νd 2
x h, the equations for the error are

   L a e 1 a = -ν∂ 2 x u in Ω 2 , e 1 a (L 2 , •) = 0, e 1 a (•, 0) = 0,    L a e 2 a = R e 1 a in Ω 2 , e 2 a (L 2 , •) = L ma e 1 a (L 2 , •), e 2 a (•, 0) = -νd 2 x h,        L ad e ad = 0 in Ω 1 , e ad (-L 1 , •) = 0, L ma e ad (0, •) = e 2 a (0, •), e ad (•, 0) = 0.
We now analyze each of the three solves separately:

First advection equation in Ω 2 : With Theorem 2.2, we find that the error e 1 a satisfies for k = 0, 1 and 2 the estimate

∂ k t e 1 a 2 L 2 x,t + |a| ∂ k t e 1 a (0, •) 2 L 2 t ≤ C ν 2 ∂ k t ∂ 2 x u 2 L 2 x,t + ∂ k t e 1 a (•, 0) 2 L 2 x . (4.5)
The case k = 0 yields the first result of the theorem. We further compute

∂ t e 1 a (•, 0) = -νd 2 x h, ∂ 2 t e 1 a (•, 0) = ν(ad x + c)d 2 x h -νd 2 x ∂ t u(•, 0), with ∂ t u(•, 0) = f (•, 0) -(ad x h + ch -νd 2 x h), so that ∂ t e 1 a (•, 0) 2 L 2 x ≤ ν 2 h H 2 x , ∂ 2 t e 1 a (•, 0) 2 ≤ ν 2 ( f (•, 0) 2 H 2 x + h 2 H 4
x ). We thus obtain for R e 1 a = (c + ∂ t ) 2 e 1 a the estimate

R e 1 a 2 L 2 x,t ≤ Cν 2 ( u 2 H 2,2 x,t + f (•, 0) 2 H 2 x + h 2 H 4 x ). (4.6)
Second advection equation in Ω 2 : Using again Theorem 2.2, we obtain the estimate We can thus conclude using (4.7).

e 2 a (0, •) 2 L 2 t ≤ C( R e 1 a 2 L 2 x,t + ν 2 h 2 H 2 x + L ma e 1 a (L 2 , •) 2 L 2 t ). To evaluate L ma e 1 a (L 2 , •), we observe that u(L 2 , •) = u a (L 2 , •), so that we have L ma e 1 a (L 2 , •) = -a∂ x e 1 a (L 2 , •) = ν∂ 2 x u(L 2 , •). Therefore, using (4.6), we get e 2 a (0, •) 2 L 2 t ≤ Cν 2 ( u 2 H 2,2 x,t + f (•, 0) 2 H 2 x + h 2 H 4 x + ∂ 2 x u(L 2 , •) 2 L 2 t ). ( 4
It remains to estimates ∂ 2 x u(L 2 , •) L 2 t and ∂ 2 x u L 2 (Ω2×(0,T )) . If the data is compactly supported, there is only one boundary layer, at x = -L 1 , and (see [START_REF] Métivier | Small viscosity and boundary layer methods[END_REF]) u(t, x) = U (t, x) + e a(x+L1)/ν U (t, 0) + O(ν), (note a < 0).

Here, U is the solution of the advection equation in Ω with data g 2 at x = L 2 . The norm of ∂ xx u, though not bounded in the entire interval Ω, is bounded in Ω 2 , since

d 2
x e a(x+L1)/ν 2 L 2 (Ω2) = |a| 3 2ν 3 (e 2aL1/ν -e 2a(L2+L1)/ν ) ∼ |a| 3 2ν 3 e 2aL1/ν , which tends to zero as ν goes to zero, because a < 0. Similarly the value at L 2 is bounded.

5. Numerical Experiments. We use a Crank-Nicolson scheme for the advectiondiffusion equation and an implicit upwind scheme for the advection equation. We discretize Ω := (-1, 1) with N = 64000 points, which leads to a spatial step ∆x = 3.125 × 10 -5 and the time step ∆t = ∆x. We choose c = 1, g 1 ≡ g 2 ≡ 0, T = 1 and the right hand side, shown in Figure 5.1 on the left, is f (x, t) = f 1 (t)f 2 (x, t), f 1 (t) = (sin 4 (4π(t -t 0 )) + sin 4 (2π(t -t 0 ))/2)χ t>t0 , t 0 = 0.1, f 2 (x, t) = e -100x 2 /4 + e -100(x-t/4-0.4) 2 + e -100(x+t/2+0.4) 2 .

Positive advection.

We choose a = 1, with the initial condition, shown in Figure 5.1 on the right, u 0 (x) = e -100(x-x0) 2 , with x 0 = -0.6. ad is very close to the viscous solution. This solution is improved with the second iteration when u a is also improved. . Right: uua L 2 x,t 5.2. Negative advection. We now consider a negative advection example, a = -1, with initial condition u 0 (x) = e -100(x-x0) 2 , with x 0 = 0.5. Figure 5.4 shows the L 2 space-time error between the viscous solution and the solution of the factorization algorithm (2.6), and also a comparison to the errors of the other coupling algorithms from the literature; the variational coupling conditions for a < 0 are -ν∂ x u ad (0, •) + au ad (0, •) = au a (0, •), and the non-variational ones are u ad (0, •) = u a (0, •). Once again the error in Ω 2 is O(ν) for each algorithm, since each algorithm solves the same advection equation in Ω 2 . However the factorization algorithm solves then a second advection equation which provides a better boundary value for the advection-diffusion problem in Ω 1 and thus can provide an error O(ν 2 ), whereas the other algorithms only give an approximation O(ν) in Ω 1 .

6. Conclusions. We introduced a new algorithm to solve advection diffusion problems with pure advection approximation in a subregion. We call this algorithm factorization algorithm, because it is based on a factorization of the underlying operator. We proved rigorous error estimates that show that our new algorithm gives solutions that are closer to the fully viscous solution of interest than other coupling algorithms in the literature. Our numerical experiments indicate that our estimates are sharp, an issue we are currently investigating using multiscale expansions.

. 11 )

 11 Let M b be the spatial part of the operator L b , i.e. L b = ∂ t + M b . We denote the boundary point where the flux enters the domain by x -, the other boundary point by x + , and define the characteristic time τ (x) :

Fig. 5 . 1 .

 51 Fig. 5.1. Left: contour plot for the right hand side in space and time. Right: initial condition a > 0

Figure 5 .

 5 Figure 5.2 shows first snapshots in time of the right hand side, and then of the viscous solution (3.1) and the solution obtained by the factorization algorithm (2.5) after one and two iterations when ν = 10 -3 . We see that in the first iteration the solution u 1ad is very close to the viscous solution. This solution is improved with the second iteration when u a is also improved.

Fig. 5 . 3 .Fig. 5 . 4 .

 5354 Fig. 5.3. Errors for a > 0 as the viscosity becomes small for our factorization algorithm compared to other coupling algorithms from the literature. Left: uu ad L 2 x,t . Right: uua L 2 x,t

) and gives a comparison to algorithms from the literature. These algorithms solve an advection-diffusion equation L ad u ad = f in Ω 1 and an advection equation L a u a = f in Ω 2 , and use for a > 0 either nonvariational transmission conditions ∂ x u ad (0, •) = ∂ x u a (0, •) and u ad (0, •) = u a (0, •), see [START_REF] Dubach | Contribution à la Résolution des Équations fluides en domaine non borné[END_REF][START_REF] Gastaldi | On the coupling of hyperbolic and parabolic systems: analytical and numerical approach[END_REF], or variational transmission conditions ν∂ x u ad (0, •) = 0 and u ad (0, •) = u a (0, •), see [START_REF] Gastaldi | On the coupling of hyperbolic and parabolic systems: analytical and numerical approach[END_REF][START_REF] Gastaldi | On the coupling of two dimensional hyperbolic and elliptic equations : Analytical and numerical approach[END_REF]. We see that the variational transmission conditions do not need an iteration in this case, one can first solve advection-diffusion, and then advection. The error is however O(ν ), as good as with one iteration of the factorization algorithm, but one needs to iterate and choosing a good relaxation parameter to ensure convergence is not easy; we chose heuristically θ = 1 450 √ ν in our computations. In the inviscid subregion Ω 2 , the error of all methods is O(ν), only the initialization step in the factorization algorithm has an error of O(1), as predicted by Theorem 3.5.