F E Kennedy

J H Currier

S Plumet

J L Duda

D P Gestwick

J P Collier

B H Currier

M-C Dubourg

Contact Fatigue Failure of Ultra-High Molecular Weight Polyethylene Bearing Components of

Introduction

Arthritis of the knee is a crippling disease which can severely limit the mobility of those afflicted by it. Joint replacement represents an appealing means to relieve the pain and increase the mobility of arthritic patients. Nearly 200,000 Total Knee Replacements (TKR) are implanted annually in the United States [START_REF] Cameron | Bone Implant Interface[END_REF]. Most TKR designs incorporate a polished metallic femoral component, usually made from a cobaltchrome alloy, in contact with a tibial bearing made from Ultra-High Molecular Weight Polyethylene (UHMWPE). An example of a typical TKR design is shown in Fig. 1. In contrast to a total hip replacement, which consists of a conformal ball-in-socket joint, the bearing in a TKR is primarily nonconformal and features a relatively concentrated contact in most designs. Whereas the kinematic contact conditions in a hip joint result in multidirectional oscillatory sliding, those in the knee are a combination of oscillatory rolling and sliding. The loading conditions in knee joints are cyclic in nature, and they vary with the flexion and extension angle. One model of the loads seen by a knee during a typical walk cycle is shown in Fig. 2 [START_REF] Iso | Test Methods for Laboratory Simulation of In-Vivo Wear of Total Knee Prostheses[END_REF]. Owing to inertia forces during the human gait, loads of more than 3-4 times body weight can be supported by the knee during a typical walk cycle (Fig. 2), and other research has shown that loads can exceed 6 times body weight in many cases [START_REF] Seireg | The Prediction of Muscular Load Sharing and Joint Forces in the Lower Extremities During Walking[END_REF]. These loads cause high contact pressures in the concentrated contact between the metallic femoral component and the UHMWPE tibial bearing (Mc-Namara et al., 1994). Therefore, after implantation the nonconforming TKR bearings are subjected to many cycles of high contact stress during oscillatory rolling/sliding contact.

Although many implanted knee prostheses have withstood Ͼ10 or even 20 years of activity without substantial surface damage (Fig. 3(a)), many other tibial bearings of TKR have shown evidence of significant wear and contact fatigue damage, requiring early replacement (Ͻ5 years after implantation) (Fig. 3(b)). The replacement procedures are costly and can have serious health consequences.

The objectives of this work have been to study the factors contributing to surface failures in UHMWPE tibial bearing components of total knee replacements.

Background

Ultra-high molecular weight polyethylene (UHMWPE) has been used as the primary bearing material in orthopaedic implants for more than 30 years. In most total knee replacements the tibial bearings are made from medical grade UHMWPE having an average molecular weight of 4 million or greater; they are either molded or machined to their final shape. Before the TKR components can be implanted in humans they must be sterilized, and until recently nearly all UHMWPE bearings used in the United States had been sterilized with gamma irradiation in air.

However, studies of bearings explanted from patients have shown dramatic variations in polyethylene wear and fatigue damage. Historically, these variations were attributed in part to patient variables, component designs, and in-vivo service histories [START_REF] Bartel | The Effect of Conformity, Thickness and Material on Stresses in Ultra-High Molecular Weight Components for Total Joint Replacement[END_REF][START_REF] Landy | Wear of Ultra High Molecular Weight Polyethylene Components of 90 Retrieved Knee Prostheses[END_REF]. Recent research has demonstrated that the damage seen in many explanted polyethylene bearings can be linked to sterilization by gamma irradiation in air and subsequent material degradation [START_REF] Sutula | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Hip[END_REF][START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF]. Over 1600 retrieved UHMWPE knee bearings, including both patellar and tibial bearings, have been examined in recent years at Dartmouth Biomedical Engineering Center [START_REF] Williams | The Impact of Sterilization Method on Wear in Knee Arthoplasty[END_REF]. Those components were provided by about 300 surgeons from North America and Europe, and they had been produced by a dozen different manufacturers. The bearings had been removed either during revision surgery or post-mortem and had been in service in-vivo for an average of about 5 years, with a range of 1 month to 20 years. Examination of the knee bearings showed that a variety of wear modes, including abrasion, burnishing and scratching, were observed on bearing surfaces, including those bearings with short in-vivo duration. The most prevalent surface failure mechanisms, however, particularly among those components that had been in service for a longer time, were delamination and contact fatigue cracking. 55% of the retrieved bearings had indications of delamination or cracking or both [START_REF] Williams | The Impact of Sterilization Method on Wear in Knee Arthoplasty[END_REF], and in nearly all cases the damage appeared to initiate subsurface. The incidence of damage was highest in bearings that had been sterilized with gamma radiation; 902 of the 1603 gammairradiated bearings showed cracking or delamination or both. The small number (n ϭ 32) of knee bearings that were known to have been sterilized by an alternative technique (ethylene oxide, or EtO) showed no evidence of contact fatigue damage, even after in-vivo durations in excess of 15 years [START_REF] Williams | The Impact of Sterilization Method on Wear in Knee Arthoplasty[END_REF]. When the specimen population was limited to gamma-irradiated components that had been in-vivo for more than four years, the percentage of tibial bearings showing either delamination or cracking increased to Ͼ65% [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF] or up to 80% [START_REF] Williams | The Impact of Sterilization Method on Wear in Knee Arthoplasty[END_REF]. The damage was so severe that about 20% of the bearings were worn through their entire thickness (Fig. 3(b)). These results are similar to those of [START_REF] Blunn | Wear in Retrieved Condylar Knee Arthroplasties[END_REF], who reported wear comparisons for 280 retrieved bearings of eight different designs. Delamination was seen in retrievals of all designs except one and in the four knee designs for which delamination frequency was reported, 65 per cent of bearings showed delamination.

About 400 of the retrieved UHMWPE bearings in the Dartmouth study were sliced through their thickness and then thin sections (200 m thick) were cut from the cross-sections using a microtome. Microtomed thin sections of retrieved tibial bearing components showed that 68% of the sections from bearings that had been gamma irradiated in air had a subsurface white band [START_REF] Williams | The Impact of Sterilization Method on Wear in Knee Arthoplasty[END_REF]. The white band was present in approximately 90% of the bearings that had cracked or delaminated in service [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF], so there was a very strong association between the presence of a subsurface white band and contact fatigue damage. Subsequent analysis of the white band using Fourier Transform Infrared Spectroscopy (FTIR) showed that it was characterized by oxidized polyethylene [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF]. The location of peak oxidation was 1-1.5 mm subsurface and was within the white band. Figure 4 shows a typical oxidation profile for specimens which had been gamma irradiated in air and shelf aged. Oxidation is measured by the maximum level of absorbance at the ketone peak (1717 cm Ϫ1 ), and that peak absorbance value is then divided by the specimen thickness (mm). The material thickness is roughly separated into three regions, a surface 'clear' layer, a subsurface oxidized ('white') layer, and a central region. The oxidation profile is symmetrical about the centerline of the material, and only the top half of the thickness is shown in Fig. 4. (Note: The layer boundaries shown on Fig. 4 are only approximate; there is no sharp demarcation between the layers.) The surface layer was clear and showed much less oxidation than the subsurface white layer; the center of the bearing thickness also showed relatively low oxidation. The oxidation takes time to occur, so the oxidation level is a function of the environment in which the bearing was irradiated and stored, and the length of time after irradiation [START_REF] Sutula | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Hip[END_REF]. Specimens which had been sterilized using EtO had significantly lower oxidation levels [START_REF] Williams | The Impact of Sterilization Method on Wear in Knee Arthoplasty[END_REF].

Mechanical testing was done on dumbbell-shaped test specimens [ASTM Type V] following ASTM standard test method D638 [START_REF] Astm | Standard Test Method for Tensile Properties of Plastics[END_REF]. The specimens were die cut from thin sections (200 m thick) that had been sliced from the UHMWPE bearings using a microtome [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF]. A series of the thin sections was tested to provide a depth profile of the mechanical properties. Tensile testing of the thin sections was done on a servohydraulic test machine to determine yield strength, ultimate tensile strength, modulus of elasticity and ductility (% elongation at failure). Poisson's ratio for the layers was not determined. Cyclic tensile tests were also done to determine the fatigue strength of the thin sections of UHMWPE. The static tensile tests showed that the material in the oxidized region ('white band') had significantly lower tensile strength and very substantially reduced ductility when compared with either unsterilized material or material from layers showing less oxidation [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF]. A summary of the test data is given in Table 1. It should be noted that there was some scatter in the modulus of elasticity results, so an alternative method for determining modulus of elasticity was also employed, based on oxidation measurements and on correlations by [START_REF] Kurtz | A Bilinear Material Model for UHMWPE Components[END_REF][START_REF] Kurtz | The Degradation Rate of UHMWPE in Total Joint Replacements[END_REF] between density changes in the oxidation zone and changes in the elastic modulus. Data for modulus of elasticity from the tensile tests (Table 1) were consistent with those calculations.

The oxidized band of material has lower tensile strength than the material above or below it, but the most striking characteristic of that layer is its very substantially reduced ductility. Cyclic tensile testing showed that specimens containing a white band suffered from early fatigue failure resulting from cracks which initiated in the oxidized white band. Specimens without a white band had fatigue lives at least 10 times longer than those with the white band [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF].

Based on the data reported above, it was concluded that the oxidation which resulted from gamma sterilization in air was responsible for reduced strength and ductility of the subsurface region of the UHMWPE bearings and that these deteriorated properties are related to the clinically-observed surface fatigue failures [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF].

Until recently, laboratory testing of materials for TKR bearings had not demonstrated the types of surface failure that had been seen in retrieved knee bearings. Most laboratory testing had concentrated on surface wear and not on contact fatigue as a surface failure mode [START_REF] Bartel | The Effect of Conformity, Thickness and Material on Stresses in Ultra-High Molecular Weight Components for Total Joint Replacement[END_REF][START_REF] Walker | Wear Testing of Materials and Surfaces for Total Knee Replacements[END_REF], despite the prevalence of delamination and cracking in retrieved bearings. In a recent study, [START_REF] Currier | In-Vitro Simulation of Contact Fatigue Damage Found in UHMWPE Components of Knee Prostheses[END_REF] found that when neversterilized or recently-sterilized UHMWPE components were tested in a knee simulator, no cracking or delamination were observed. However, when the bearings had been gamma sterilized in air and aged prior to testing, cracking and delamination did occur under the same conditions and it was similar to the surface damage found in retrieved knee bearings [START_REF] Currier | In-Vitro Simulation of Contact Fatigue Damage Found in UHMWPE Components of Knee Prostheses[END_REF]. Thus it was concluded that in order to duplicate the clinical damage mechanisms, it is necessary to use bearing material which had been subjected to the same sterilization and aging conditions as the material of clinical knee components.

The goal of this study was to complete the picture of contact failure of UHMWPE bearings by studying the tribological conditions responsible for the initiation and propagation of damage in gamma irradiated and aged components. The study consisted of two thrusts:

• Tribotesting of UHMWPE bearing material in a knee simulator and a rolling/sliding contact tester • Numerical analysis of the contact stress distribution in UH-MWPE tibial bearings under service conditions.

Experimental Methods and Materials

The current study employed two different tribotest machines to assess the behavior of UHMWPE bearing material. Although some investigators have used sliding wear tests to study the wear of total knee replacement bearing materials (e.g., [START_REF] Walker | Wear Testing of Materials and Surfaces for Total Knee Replacements[END_REF], those tests do not provide the same contact kinematics as in the knee, so they usually result in surface failure modes that are different from the delamination and contact fatigue observed in clinical retrievals. In this study, two different testers, a single station knee simulator and a rolling/sliding wear tester, were used to replicate the contact fatigue failure observed in retrieved tibial bearings. This study focused on identifying these fatigue-related surface failure modes, so no quantified wear measurements (e.g., mass loss) were recorded during any of the tests. The lubricant used in all tests was distilled water. Although bovine serum has physical properties that are more similar to human synovial fluid, distilled water has been found to give a reasonable simulation of contact conditions which influence subsurface-originated rolling/ sliding contact fatigue phenomena [START_REF] Currier | In-Vitro Simulation of Contact Fatigue Damage Found in UHMWPE Components of Knee Prostheses[END_REF]. Knee Simulator. A single station knee simulator was used to test the failure modes in the prosthesis. The knee simulator is described in detail elsewhere [START_REF] Currier | In-Vitro Simulation of Contact Fatigue Damage Found in UHMWPE Components of Knee Prostheses[END_REF]. The simulator uses actual TKR components, including a polyethylene tibial insert and a matching metallic femoral component. Load, flexion angle and cycle rate are the variables that can be controlled. All tests reported here were run with a flexion range of 0 to 35°. Several different load cycles were used in the test program. Although the load cycles didn't duplicate exactly the one shown in Fig. 2, the test cycles all provided a relationship between maximum load and flexion angle similar to that shown in the figure. In the most severe load cycle tested the maximum load was 2000 N (450 lb), whilst in the most gentle cycle the maximum load was 780 N (175 lb). Lubrication was provided by a bath of distilled water. A cooling loop placed in the bath kept the temperature in the range of 22-24°C.

All test specimens were tibial bearings manufactured for implantation and taken from shelf inventory for testing. All were machined from UHMWPE bar stock that had been ram extruded using GUR-415 resin powder as a starting material. The components were gamma-irradiation sterilized in an ambient air environment according to the industry standard protocol of 2.5-4.0 megarads from a Co 60 source. The bearings were shelf aged in air for periods ranging from 0 to Ͼ8 years after gamma sterilization. Tibial bearing thickness ranged from 8 -18 mm and all were of the same design.

A single matching metallic femoral component was used for all tests in this series. It was made of cast cobalt chrome (Haynes 21 alloy), and polished according to the normal protocol by the implant manufacturer. The surface of the metallic component was cleaned with a solvent cleaner before each test, but no attempt was made to remove all adhered transfer films.

Rolling/Sliding Tester. The rolling/sliding tester enabled tests to be run at a well-controlled set of constant contact conditions, in contrast to the varying load and contact kinematics of the knee simulator. The rolling/sliding wear tester employs two counter-rotating disks to simulate the rolling and sliding in tibiofemoral loading. Details of the tester are given elsewhere [START_REF] Kennedy | Failure of UHMWPE polymer surfaces of knee prostheses: a rolling sliding contact[END_REF]. All disks were 62.5 mm diameter and 25 mm thick. The polymer test disks were machined from UHMWPE bar stock that had been ram extruded using GUR-415 polyethylene resin powder as a starting material. The UHMWPE disks were sterilized by gamma irradiation in air and artificially aged in an attempt to simulate extended time on the shelf. A number of artificial aging techniques have been developed in recent years to achieve in a short period the type of long-term oxidation damage that has been observed in shelf-aged or explanted bearings (e.g., [START_REF] Saum | Oxidation Versus Depth and Time for Polyethylene Gamma Sterilized in Air[END_REF][START_REF] Sun | Development of an Accelerated Aging Method for Evaluation of Long-Term Irradiation Effects on Ultrahigh Molecular Weight Polyethylent Implants[END_REF]. Two different accelerated aging protocols were used in the study: one atmosphere of oxygen at 80°C for 17 days, and three atmospheres of oxygen at 70°C for 17 days. Control specimens were UHMWPE disks that had been neither sterilized nor aged, although some specimens which had been sterilized but not aged were also tested.

The counterface disk was machined from cobalt chrome alloy and polished to a finish similar to that of the surface of the manufactured femoral component.

The rolling/sliding tests were done at 2 Hz and an amplitude of 15 mm, with a specific slide:roll ratio of 1:5. Tests were done at normal loads ranging from 135 N to 2000 N (450 lb). A test was continued until surface damage was first observed, or to a maximum of 2 million cycles. All tests were conducted in a distilled water bath at a temperature of 24°C.

Post-test analysis of samples from both types of tests included optical microscopy, both before and after sectioning the component to a 200 m thickness with a microtome, and Scanning Electron Microscopy (SEM). FTIR examination of thin sections was also carried out for some specimens.

Experimental Results

Knee Simulator. Some of the knee simulator results were reported by [START_REF] Currier | In-Vitro Simulation of Contact Fatigue Damage Found in UHMWPE Components of Knee Prostheses[END_REF] and will be summarized here for completeness. Other test data have confirmed those findings.

Tests of unsterilized or recently-sterilized (but unaged) control specimens were run on the knee simulator to a maximum of 2.5 million cycles with a load cycle having a peak load of 2000 N. No cracking or delamination was found and the only wear that appeared was slight burnishing at the surface of the UHMWPE components. Measurements of the wear magnitude were not made, but the wear patterns were very similar to those found on retrieved bearings that did not suffer from cracking or delamination. Surfaces that had been aged long enough to develop an oxidized white layer, on the other hand, showed more severe surface and subsurface damage.

The sterilized and shelf aged components were run for significantly fewer cycles than the unsterilized or recently-sterilized controls (100,000 to 1 million, as opposed to 2.5 million cycles), yet the aged bearings exhibited more severe surface damage. All aged samples showed surface wear, including both flakes of polyethylene peeling from the surface and deep wear grooves on their articulating surfaces. The irradiated and aged specimens also showed surface delamination of the type found in retrieved tibial bearings. Delamination in the articular area appeared as white blister-like regions visible beneath the intact surface layer of the component [START_REF] Currier | In-Vitro Simulation of Contact Fatigue Damage Found in UHMWPE Components of Knee Prostheses[END_REF].

After the tests the UHMWPE specimens were sliced normal to the worn surface and thin cross-sections were made by either microtoming or slicing with a diamond saw so that the subsurface region could be examined. Examination of those cross-sections revealed that most of the aged specimens had a subsurface oxidized band. The severity of the oxidation increased with the increasing aging duration after gamma irradiation. The amounts of delamination and subsurface cracking were found to increase as the severity of the subsurface oxidation increased (Table 2). Cracks in the oxidized band were more numerous beneath the contact zone on the bearings than in the non-articular region; cyclic loading and unloading appeared to affect the size and severity of the defects and cracks [START_REF] Currier | In-Vitro Simulation of Contact Fatigue Damage Found in UHMWPE Components of Knee Prostheses[END_REF]. Subsurface cracks in the components with the most severe oxidation band were quite severe and had propagated nearly to the surface (Fig. 5). The subsurface cracks seen in Fig. 5 are very similar to those observed in explanted prostheses (e.g., [START_REF] Blunn | The Dominance of Cyclic Sliding in Producing Wear in Total Knee Prostheses[END_REF].

Rolling/Sliding Tester. The polymer specimens tested on the rolling and sliding tester produced varying degrees of wear characteristics, including burnishing, surface cracking, surface pitting, subsurface cracking and delamination. It was found that the amount of damage and the depth at which the damage developed were related to the irradiation and aging methodology (Table 3). In general, specimens that had not been sterilized showed little surface damage, even after up to 2 million cycles of oscillation. The same was true of specimens which had been sterilized but not aged after sterilization. For both types of specimens there was some light surface wear or burnishing and a small number of surface originated cracks were found on specimens which had been tested for at least one million cycles. When cracks were found they were oriented perpendicular to the sliding direction and extended less than 1 mm beneath the surface. Those results were in agreement with earlier rolling/sliding tests of unsterilized UHMWPE disks [START_REF] Kennedy | Failure of UHMWPE polymer surfaces of knee prostheses: a rolling sliding contact[END_REF].

Tests of the aged specimens (SA1 or SA3) generally resulted in cracks and/or pits on the surface, as long as the test had been carried out for a sufficient number of cycles or at high enough load (Table 3). The facts that cracks and pits did not appear in tests with low normal load and that at a given load there were more pits and cracks after a greater number of cycles are both consistent with a fatigue-type phenomenon. This is confirmation that the surface failure of the sterilized and aged UHMWPE components is related to contact fatigue.

Specimens which had been aged in one atmosphere of oxygen after sterilization (SA1) were found to have cracks and pits that originated at or very close to the contact surface. An example is shown in Fig. 6. On the other hand, specimens that were aged in three atmospheres of oxygen (SA3) showed evidence of subsurface-initiated cracks and spalls. The subsurface cracks were best viewed by slicing the specimens parallel to the sliding direction and viewing the subsurface region in cross-section. One such example is shown in Fig. 7. It is evident that a crack had initiated subsurface and had propagated approximately parallel to the surface at a depth of 300m. Reasons for the differences between the SA1 and SA3 specimens became clear after they were analyzed for oxidation. Microtomed thin specimens of the disk specimens were analyzed using Fourier Transform Infrared Spectroscopy (FTIR) according to procedures that had been developed earlier for analyzing tibial bearings [START_REF] Sutula | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Hip[END_REF]. Oxidation profiles (ketone levels) for SA1 and SA3 specimens are shown in Figs. 8 and9, respectively. It is apparent from Fig. 8 that oxidation in the SA1 specimens reached a peak at the surface and decreased rapidly in the subsurface region. On the other hand, the SA3 specimens had an oxidation peak beneath the surface (Fig. 9), but there was more surface oxidation than in the shelf-aged tibial bearings (Fig. 4). From these oxidation profiles, it can be concluded that in all cases the contact fatigue cracks and delamination began in the most heavily oxidized region. This region was located at a depth of about 1-2 mm for the shelf-aged bearings tested in the knee simulator, whereas it was at a depth of about 500 m for the SA3 specimens and it was at the surface for the SA1 specimens. Although there was less surface oxidation in the SA3 specimens than in the SA1 specimens, there was still considerably more surface oxidation in all accelerated-aged specimens than in shelf-aged bearings. Surfaceoriginated damage was most prevalent with the SA1 specimens, and least prevalent with the shelf-aged knee simulator specimens; therefore surface-originated damage tended to occur when there was surface oxidation.

Contact Stress Analysis

Methodology. In order to understand the conditions in the region beneath the contact where the cracks initiated, a contact stress analysis was carried out. Other investigators have studied the stress distribution in UHMWPE tibial bearings in recent years, and most have used finite element methods. For example, Sathasivam and Walker (1998) did a three-dimensional elastic finite element analysis to determine contact area, contact pressure distribution and subsurface strain energy density accumulation to predict the probability of subsurface damage in tibial bearings. [START_REF] Estupinan | Residual Stresses in Ultra-High Molecular Weight Polyethylene Loaded Cyclically by a Rigid Moving Indenter in Nonconforming Geometries[END_REF] and [START_REF] Reeves | A Two-Dimensional Model of Cyclic Strain Accumulation in Ultra-High Molecular Weight Polyethylene Knee Replacements[END_REF] did two-dimensional elastoplastic analyses of UHMWPE bearings in order to study the residual stresses and plastic strain accumulation in the contact region of the bearings. None of these analysis, however, have accounted for the fact that gamma irradiated and aged tibial bearings are inhomogenous as a result of the changes in mechanical properties that occur with oxidation. In fact, finite element methods generally require large computer times to achieve accurate stress distributions within three-dimensional non-conformal contacts with layered material properties. To overcome that problem, a combined analytical/numerical method was developed by [START_REF] Plumet | Modelisation d'un milieu multicouche 3D sous sollicitations de contact: Application aux prosthe `ses de genou ste ´rilise ´es[END_REF] based on integral transform techniques. That method was used in this work to study stresses in knee prosthesis bearings with layered material properties resulting from oxidation.

Details of the analytical/numerical method used here are presented elsewhere [START_REF] Plumet | Modelisation d'un milieu multicouche 3D sous sollicitations de contact: Application aux prosthe `ses de genou ste ´rilise ´es[END_REF][START_REF] Plumet | A 3D Model for a Multilayered Body Loaded Normally and Tangentially against a Rigid Body: Application to Specific Coatings[END_REF]. The method is composed of the following steps:

-The UHMWPE bearing is assumed to be composed of n layers. Each layer is assumed to be isotropic, linearly elastic, smooth and homogeneous, and is characterized by its elastic constants (E i , i ) and thickness e i . The proper- ties in each layer were those determined in the tensile tests (Table 1). -The contacting femoral condyle is assumed to be rigid and is characterized by two radii of curvature (Fig. 10). The radii of curvature are a function of the flexion angle (Fig. 1). -The model is based on linear elasticity theory, unilateral contact with friction, and Fourier transforms. Fast Fourier Transforms were used to enable computational time to be reduced. -Boundary conditions are applied as shown in Figure 10.

The bottom surface was assumed to be fixed to a 'rigid' metallic backing. Contact between the condyle and the bearing was assumed to be frictionless. -An iterative algorithm developed by [START_REF] Kalkar | Three Dimensional Elastic Bodies in Rolling Contact[END_REF] is used to determine contact area, sliding conditions and tractions within that area, and contact pressure distribution between the rigid condyle and the tibial bearing surface. -The surface tractions are used to determine the stress field within the multilayered body. -The analysis was carried out at four different positions (flexion angles) and the contact radii and normal loads consistent with those positions were used.

Results.

A large number of analyses were carried out to determine the influence of flexion angle, normal load, thickness of tibial bearing, interfacial friction and layer properties on the contact area, contact pressure and stress distribution. Those results are given elsewhere [START_REF] Plumet | Modelisation d'un milieu multicouche 3D sous sollicitations de contact: Application aux prosthe `ses de genou ste ´rilise ´es[END_REF]; only the ones most relevant to the oxidized tibial bearing are presented here.

Analyses of homogeneous tibial bearings showed that, although the maximum normal stress was the contact pressure on the surface, the largest von Mises equivalent stress and maximum shear stress occurred within the bearing beneath the center of the contact, consistent with Hertzian contact stress analysis [START_REF] Johnson | Contact Mechanics[END_REF]. Properties for the homogeneous model were those for unsterilized UHMWPE (Table 1).

When an oxidized layer with higher modulus was included in the model, it led to smaller contact area, larger contact pressure, and an increase in von Mises equivalent stress (Table 4). For the 5-layer model, the five layers included surface 'clear' layer, subsurface oxidized layer, center region, subsurface oxidized layer, and surface 'clear' layer. Properties for each of those layers were those given in Table 1. Each layer was assumed to have the same Poisson's ratio ( ϭ 0.3). Other conditions are given in Table 4.

The distribution of von Mises equivalent stress in a five-layer model of a tibial bearing is shown in Fig. 11 for a given normal load (1000 N), bearing thickness (8.0 mm), radii of curvature (Rx 1 ϭ 60 mm, Rx 2 ϭ 130 mm), and for a semi-bent flexural position of the knee. [It should be noted that both the prosthesis geometry (radii of curvature) and normal load are different for the case in Fig. 11 than for the one analyzed in Table 4.] It is apparent that the maximum von Mises stress occurs in the oxidized layer, which is also the layer with the lowest strength and by far the lowest ductility. This is the most probable reason for the initiation and propagation of the subsurface cracks found in retrieved tibial bearings and bearings tested in the knee simulator (Fig. 5).

It was found [START_REF] Plumet | Modelisation d'un milieu multicouche 3D sous sollicitations de contact: Application aux prosthe `ses de genou ste ´rilise ´es[END_REF] that the maximum von Mises stress increases with larger normal load, with greater knee flexion (corresponding to smaller radius of the contacting tibial component, Fig. 1), and with a decrease in the depth of the stiffer oxidized layer. It was noted above that the oxidized layer was on the surface for accelerated aging protocol SA1, was slightly subsurface (about 0.5 mm) for specimens aged by protocol SA3, and was deepest (1-2 mm) for shelf aged or in-vivo aged bearings. Therefore, for given normal load and component radii, the maximum von Mises stress (or maximum shear stress) would be greater for the SA1 specimens than for SA3, and both accelerated aging cases would produce a larger von Mises stress than for bearings that had been shelf aged after gamma irradiation. Therefore, fatigue lives would be expected to be lower for the accelerated aged specimens than for shelf aged material under similar contact conditions.

The property of the oxidized material which has the most influence on the stress distribution in the layered model is its modulus of elasticity. Modulus information was not gathered for the aged specimens, so the depth of maximum shear stress shown in Fig. 7 is only approximate. However, it is clear that the von Mises stress (or maximum shear stress) will be quite high in the oxidized band, whether or not the peak value of those stresses occurs within the band (Fig. 11). It is the relatively large cyclic stress in the region of much reduced ductility and strength that is the cause of contact fatigue cracks in the oxidized band.

Another characteristic of the stress distribution in the oxidized tibial bearing which could be important is the interfacial shear stress between the layers. That interfacial shear stress could influence crack propagation at the interface, leading to delamination. It was found [START_REF] Plumet | Modelisation d'un milieu multicouche 3D sous sollicitations de contact: Application aux prosthe `ses de genou ste ´rilise ´es[END_REF] that the greatest interfacial shear stress occurs at the interface between oxidized layer and the clear layer above it, leading to the possibility of delamination at those interfaces. That result is also in accord with observations of delamination in retrieved bearings and those tested in the knee simulator.

Conclusions

Tibial bearings of total knee prostheses made from ultrahigh molecular weight polyethylene (UHMWPE) which has been gamma irradiated in air are susceptible to oxidation when aged in oxygen-containing environments. The oxidation reaches a peak about 1 to 2 mm beneath the surface when aged on the shelf or in-vivo. The peak oxidation takes four years or more to develop (in air or in-vivo). Accelerated aging protocols have been developed which enable oxidation to develop in a much shorter time, but the depth of the peak oxidation may differ from that in shelf-aged components, and the peak oxidation could occur on the surface under some conditions. The material properties in the oxidized layer are significantly altered, including a very substantial reduction in ductility, reductions in fatigue strength and tensile strength, and modest increases in density and modulus of elasticity.

The stress distribution beneath the contact region in a tibial bearing is affected by the presence of the oxidized layer, and the largest von Mises stress (or maximum shear stress) is generally higher than in unoxidized material. That peak stress tends to occur within the oxidized layer in naturally aged material. In accelerated aged material the peak stress may occur within or beneath the oxidized layer, depending on the location of the peak oxidation.

The high stress and reduced ductility and strength in the oxidized layer lead to the initiation of fatigue cracks in that layer. This can result in delamination, spalling, and contact fatigue failure of the UHMWPE bearings. Contact fatigue failure of UHMWPE bearings has been found in gamma irradiated and aged tibial bearings in a knee simulator. Those failures were very similar to those which have been noted in explanted tibial bearings. Similar contact fatigue failures have also been noted in oxidized UHMWPE material which has been tested in a rolling/sliding contact tester. The failures in material which has been aged in an accelerated manner are similar in nature to those which occurred in shelf aged bearings in the body or in the knee simulator. The difference was the shallower depth of the crack initiation in the accelerated aged material than in shelf aged bearings, and the reason for that was the fact that the oxidized layer was closer to the surface after accelerated aging than after natural (shelf) aging. No matter what the depth of the oxidized layer, the fatigue cracks appeared to always initiate in that layer.

These results strengthen the conclusions of the earlier work by [START_REF] Collier | Impact of Gamma Sterilization on Clinical Performance of Polyethylene in the Knee[END_REF]. The clinical implications are so important that all American manufacturers of total knee replacements have recently stopped using gamma irradiation in air to sterilize their UHMWPE bearings.

Fig. 3

 3 Fig. 3(b) Retrieved tibial bearing which experienced catastrophic surface damage in the form of cracking and delamination within 8 years of service

Fig. 4 FTIR

 4 Fig. 4 FTIR-measured oxidation level (ketone peak absorbance level/ thickness in mm) versus depth (microns) for a never-implanted polyethylene tibial component 8 years after gamma sterilization in air

Fig. 5

 5 Fig. 5 Thin cross-section of tibial component tested in knee simulator for 150,000 cycles at a maximum normal load of 1225 N (275 lb.) in distilled water. It is clear that fatigue cracks have initiated at a depth of about 2 mm subsurface but have not yet propagated to the surface.

Fig. 8 FTIR

 8 Fig. 8 FTIR-measured oxidation level (ketone peak level/thickness in millimeters) versus depth (millimeters) for UHMWPE specimens that had been aged in one atmosphere of oxygen at 80°C for 17 days after gamma sterilization in air (SA1)

Fig. 11

 11 Fig. 11 Distribution of von Mises equivalent stress beneath contact with five-layer model of a particular sterilized and oxidized tibial bearing for a given normal load (1000 N), bearing thickness (8.0 mm), radii of curvature (Rx 1 ؍ 60 mm, Rx 2 ؍ 130 mm), and for a semi-bent flexural position of the knee

Table 1 Results of tensile tests of thin sections from retrieved UHMWPE bearings (adapted from Collier et al., 1996 and Currier et al., 1998)

 1

		Gamma sterilized in air and aged
		Near-surface	Subsurface		Unsterilized
		'clear' layer	'white' band Center	(benchmark)
	Ultimate tensile
	strength
	(MPa)	35	20	27	55-60
	% Elongation
	at fracture	230	15	70	550
	Modulus of
	elasticity
	(MPa)	600	1000	800	800

Table 2 Effect of oxidation level on results of knee simulator tests. All tests were run in distilled water at 22°-24°C with maximum load of 1225 N.

 2

	Subsurface
	oxidation	Test duration		Subsurface
	level	(cycles)	Surface wear	cracking
	None	500,000	Wear grooves and	None
			surface flaking
	Moderate	500,000	Wear grooves and	Many small cracks
			surface flaking	in oxidized
				layer
	Severe	150,000	Mild wear	Severe cracking
			grooves and
			surface flaking

Table 3 Results of the rolling-sliding tests. All tests were run in distilled water at 22°-24°C.

 3

		Normal
	Sterilization	load	Test duration		Subsurface
	condition	(N)	(cycles)	Surface wear	cracking
	Unsterilized	900 N	1 million	Burnishing	No
	Unsterilized 2000 N	2 million	Burnishing, some	No
				surface cracks
	S	900 N	1 million	Minor surface wear	No
	SA1	135 N	140,000	Slight burnishing	No
	SA1	630 N	140,000	Surface cracks,	No
				small pits
	SA1	630 N	65,000	Some surface	No
				cracks
	SA1	690 N	930,000	Many surface	No
				cracks, deep pits
	SA1	1040 N	140,000	Surface cracks,	No
				some pits
	SA3	900 N	14,000	Subsurface cracks,	Yes
				few spalls
	SA3	900 N	22,000	Subsurface cracks,
				spalls

Yes

S ϭ Gamma sterilized in air, not aged. SA1 ϭ Gamma sterilized in air, aged in of one atmosphere of oxygen at 80°C for 17 days SA3 ϭ Gamma sterilized in air, aged in three atmospheres of oxygen at 70°C for 17 days

Table 4 Influence of subsurface oxidation on contact pres- sure and on the magnitude and depth of maximum von Mises stress

 4 A single prosthesis geometry and knee flexion position is considered, with applied normal load ϭ 1500 N, total bearing thickness ϭ 8 mm, and radii of curvature ratio Rx 1 /Rx 2 ϭ 0.18

		Max. contact	Max. von Mises
		pressure	stress vm max	Depth of vm max
	Model	(MPa)	(MPa)	(mm)
	Unsterilized	40.3	12.6	1.6
	Sterilized
	(five-layer)	40.5	15.0	1.8

Acknowledgment

The authors gratefully acknowledge the support of DePuy Incorporated, a Johnson and Johnson Company, and the Veterans Administration Research and Rehabilitation Service through grant #A-473-3DA.


Konstantin Brenner

Acceleration of Newton's method using nonlinear Jacobi preconditioning

Keywords: Mildly nonlinear systems, Newton's method, Jacobi-Newton method, nonlinear preconditioning, monotone convergence MSC (2010): 58C15, 65H10, 65H20, 65M22

For mildly nonlinear systems, involving concave or convex diagonal nonlinearities, semiglobal monotone convergence of Newton's method is guaranteed provided that the Jacobian of the system has a nonnegative inverse. However, regardless of this convergence result, the efficiency of Newton's method becomes poor for stiff nonlinearities. We propose a nonlinear preconditioning procedure inspired by the Jacobi method and resulting in a new system of equations, which can be solved by Newton's method much more efficiently. The obtained preconditioned method is shown to by globally convergent.

Introduction

Let N be a positive integer we consider the problem of finding u ∈ R N satisfying

f (u) + Au = b, (1)
where A belongs to the set of real N × N matrices, denoted in the following by M(N ), and f is a diagonal mapping given by

f : u →    f 1 (u 1 ) . . . f N (u N )    .
Because of the applications that we have in mind we will assume that f i are only defined on R ≥0 . More specifically, the analysis presented in this article will be based on the following assumptions:

(A 1 ) For each 0 ≤ i ≤ N , the function f i is a continuous bijection from R ≥0 to R ≥0 belonging to C 1 (0, +∞). The matrix A has zero diagonal and nonpositive off-diagonal elements, and for any u > 0 the inverse of f (u) + A exists and is nonnegative. We assume in addition that b ∈ R N ≥0 .

(A 2 ) For each 0 ≤ i ≤ N , the function f i is concave.

Let us remark that the assumption (A 1 ) implies that f (u) + A is an M-matrix and therefore has a positive diagonal (see e.g. 2.4.8 of [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF]); it follows in turn that f is increasing and therefore f (0) = 0. We also remark that the derivatives of f i are potentially unbounded at the origin; we will denote f i (0) = lim u→0+ f i (u). Finally, we remark that the analysis presented in the article can be trivially adapted to the case of f i being convex instead of concave.

The system (1) can be found in the numerical modeling of flow and transport processes. In particular it arises from the discretization of the nonlinear evolutionary PDEs of the form

∂β(u) + div (Vu -λ∇u) = γ(u), (2)
where V is some given velocity field. Applying the backward Euler scheme and some space discretization method to (2) one typically gets the discrete problem of the form

β(u n h ) -β(u n-1 h ) ∆t + M -1 Su n h = γ(u n h ) + σ n h , (3)
where u n h , u n-1 h are the vectors of the discrete unknowns associated with two sequential time steps, while M and S are respectively the mass and the stiffness matrices, and the vector σ n h represent the effect of the boundary conditions.

To fix the ideas let's assume that the Dirichlet boundary conditions are imposed. Several space discretization methods ensure (possibly under some geometrical condition on the mesh) that the matrix M -1 S is an M-matrix. In the presence of diffusion (that is λ > 0), the examples of such monotone discretization schemes would be the standard finite volume method with two-point flux approximation, or P 1 finite element method with mass lumping under the Delaunay condition on the underlying mesh (see [START_REF] Eymard | Finite Volume Methods, Handbook of Numerical Analysis[END_REF]). Let us mention that the monotone discretizations are not only beneficial to the nonlinear solver (as it is going to be discussed in this paper), but also allow to preserve the local maximum principle on the discrete level, thus avoiding any spurious oscillations of the discrete solution.

Let D denote the diagonal of M -1 S and let A = ∆t M -1 S -D . Setting

f (u) = β(u) + ∆t(Du -γ(u))
the system (3) can be written in the form of (1) as

f (u n h ) + Au n h = β(u n-1 h ) + ∆tσ n h .
Given the assumptions (A 1 ) -(A 2 ) on the mapping f , and thus on the nonlinearities β(u) and γ(u), several physical models are relevant. Such models are for example the porous media equation [START_REF] Vázquez | The Porous Medium Equation -Mathematical theory[END_REF], models of transport in porous media with adsorption (using e.g. the Freundlich isotherm [START_REF] Bear | Modeling groundwater flow and pollution[END_REF]), the Richards' equation [START_REF] Duijn | Nonstationary filtration in partially saturated porous media[END_REF], [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF], or the Dupuit-Forchheimer equation [START_REF] Bear | Modeling groundwater flow and pollution[END_REF] (provided that convection is discretized using an explicit scheme). Let us further remark that the analysis and the algorithms presented in this paper can be extended to the Hele-Shaw or Stefan like problems (see Remark 2.3 below), where β(u) is no longer a function, but rather a monotone graph of the form

β(u) = ζH(u) + β(u),
where β is a nondecreasing C 1 concave function, ζ is a nonnegative real number and H denotes the multivalued Heviside graph. In [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF] this type of nonlinearity has been addressed through the parametrization of β, that is a couple of the functions τ → (u(τ ), v(τ )) with v(τ ) ∈ β(u(τ )) for all τ . Then, the problem has been reformulated in terms of this new variable τ .

Due to its quadratic convergence in the vicinity of a solution, Newton's method is a very popular tool for solving systems of algebraic equations, and in particular, those arising from the discretization of the nonlinear PDEs. Let F be some mapping from R N to R N and assume that F is differentiable in some appropriate sense. Then, starting with some initial guess u 0 ∈ R N , Newton's method generates the sequence (u n ) n defined by

u n+1 = u n -F (u n ) -1 F (u n ), n ≥ 0, (4)
which hopefully converges to some u s ∈ R N satisfying F (u s ) = 0. Unfortunately the sequence (u n ) n may not converge, in particular the celebrated Newton-Kantorovich theorem [START_REF] Kantorovich | On Newtons method for functional equations[END_REF], [START_REF] Ortega | The Newton-Kantorovich theorem[END_REF] ensures convergence of (u n ) n only if the initial guess u 0 is sufficiently close to u s . To overcome this limitation, multiple modifications of a basic Newton's method, involving line search, trust region or homotopy continuation, have been proposed [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF], [START_REF] Jr | Numerical methods for unconstrained optimization and nonlinear equations[END_REF].

For system (1) and under Assumption (A 1 ) there are some variants of Newton's method that ensure convergence of (u n ) n under some mild if any assumptions on u 0 . Those algorithm would typically generate a monotone sequence of lower or upper solutions converging to u s . Let us briefly review some of those monotone methods. First of all we remark that under assumption (A 2 ) the sequence generated by (4) will converge monotonically toward any positive solution u s as soon as the initial guess u 0 satisfies F (u 0 ) ≤ 0 (see Proposition 2.3 below), in particular the sequence (u n ) n satisfies u n ≤ u n+1 ≤ u s for all n ≥ 0. This semi-global convergence result follows from a more general Monotone Newton Theorem (MNT), which were originally introduced by Baluev [START_REF] Baluev | On the method of Chaplygin[END_REF] and is derived from two major assumptions: F is either convex or concave, and F (u) -1 is nonnegative. In this article we use a slightly weaker version of this theorem (Theorem 2.1 below), for more general results we refer to [START_REF] Ortega | Monotone iteration for nonlinear equations with applications to Gauss-Seidel methods[END_REF], [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF], [START_REF] Potra | On the monotone convergence of Newton's method[END_REF] and [START_REF] Potra | Newton-like methods with monotone convergence for solving nonlinear operator equations[END_REF].

If the left hand side of (1) is neither convex nor concave, then the original Newton's method can be modified in a way that the monotone convergence is preserved. This can be achieved either by employing a so-called method of the accelerated monotone iterations [START_REF] Pao | Accelerated monotone iterative methods for finite difference equations of reactiondiffusion[END_REF], [START_REF] Pao | Accelerated monotone iterations for numerical solutions of nonlinear elliptic boundary value problems[END_REF], or by means of a nested Newton's method [START_REF] Brugnano | Iterative Solution of Piecewise Linear Systems and Applications to Flows in Porous Media[END_REF], [START_REF] Casulli | A Nested Newton-Type Algorithm for Finite Volume Methods Solving Richards' Equation in Mixed Form[END_REF], [START_REF] Casulli | Iterative solutions of mildly nonlinear systems[END_REF].

The accelerated monotone iterations, presented in [START_REF] Pao | Accelerated monotone iterations for numerical solutions of nonlinear elliptic boundary value problems[END_REF], make use of both lower and upper solutions of F (u) = 0. In particular, given u 0 and u 0 that satisfy F (u 0 ) ≤ 0 ≤ F (u 0 ), it generates a couple of sequences (u n ) n and (u n ) n defined by max

u n ≤ξ≤un F (ξ) (v n+1 -v n ) + F (v n ) = 0 n ≥ 0, (5)
with v standing either for u or u. One shows that u n (resp. u n ) is lower (resp. upper) solution of F (u) = 0, that both sequences satisfy u n ≤ u s ≤ u n for all n ≥ 0 and converge monotonically toward u s . The latter inequality obviously provides a useful error estimate. Note that, for each n ≥ 0, one has to solve two linear systems resulting from ( 5) with v = u and v = u. However, those systems only differ by their right-hand-sides, and this situation can be efficiently handled by some linear solvers.

The second method originated from [START_REF] Brugnano | Iterative Solution of Piecewise Linear Systems and Applications to Flows in Porous Media[END_REF], is based on some particular splitting F (u) = F 1 (u) -F 2 (u), where both mappings F 1 and F 2 are either concave or convex. The system F 1 (u) -F 2 (u) = 0 is solved trough a nested iterative linearization process. The outer loop of the method generates the sequence (u n ) n defined through the following partial linearization scheme

F 1 (u n )(u n+1 -u n ) + F 1 (u n ) -F 2 (u n+1 ) = 0, n ≥ 0. (6)
One shows that the solution to [START_REF] Eymard | Finite Volume Methods, Handbook of Numerical Analysis[END_REF] exists and that the sequence (u n ) n monotonically converges to u s . For each n ≥ 0, the nonlinear system (6) can be solved again by Newton's method. This results in an inner loop generating a sequence that monotonically converges to u n+1 . Remark that each inner iteration of the algorithm requires solving a linear system, therefore the total count of linear solves is N outer × N inner .

In contrast with the aforementioned methods we do not aim to relax the convexity/concavity assumption in MNT. Instead, our objective is to accelerate the convergence of the algorithm (4). This will be achieved through a nonlinear preconditioning procedure that preserves the structure required by MNT. As a side note, we remark that, in principle, the preconditioning proposed in this article can be combined with the modified Newton's methods from [START_REF] Pao | Accelerated monotone iterations for numerical solutions of nonlinear elliptic boundary value problems[END_REF] and [START_REF] Casulli | Iterative solutions of mildly nonlinear systems[END_REF].

To motivate our study, let us remark that, despite the monotone convergence result, the efficiency of Newton's method applied to (1) can be very poor especially for stiff problems with f (0) = +∞. To give an example let γ(u) = 0 and β(u) = u m , m > 1 (this choice corresponds to the porous media equation [START_REF] Vázquez | The Porous Medium Equation -Mathematical theory[END_REF]), we demonstrate in the numerical section 3 that the convergence of Newton's method can be very slow; moreover the required number of iterations increases with m. The numerical experiment also demonstrates that the efficiency of Newton's method can be greatly improved by a simple change of the variable u = β(v). Let us note that for Richards-like parabolic-elliptic problems with β (u) = 0 for u ≥ u s > 0 the similar change-of-variable trick can be performed using the variable switching technique as suggested in [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF]. Compared to the initial formulation of (1) the drawback of the change-of-variable approaches is that the concavity of the problem is lost, and therefore the monotonic convergence is no longer guaranteed.

In this article, we reformulate (1) in a way that preserves the concavity of the system while offering a much faster convergence of the nonlinear solver. Since the modified system is similar to one obtained in the Jacobi method, we refer to our approach as the Jacobi-Newton method, or the Jacobi preconditioned Newton's method. Note that the Jacobi method can be viewed as a domain decomposition method with the minimal subdomain size and the minimal algebraic overlap. In this regard our approach can be related to the nonlinear domain decomposition methods (see e.g. [START_REF] Cai | Nonlinearly preconditioned inexact Newton algorithms[END_REF] and [START_REF] Dolean | Nonlinear preconditioning: how to use a nonlinear Schwarz method to precondition Newton's method[END_REF]).

Because the mapping f is diagonal, strictly increasing and continuous it admits an inverse for all u ≥ 0. Let g be a diagonal mapping which coincides with f -1 on R N ≥0 , we consider the following left and right-preconditioned problems

F l (u) := u -g(b -Au) = 0 (7) and F r (ξ) := ξ + Ag(ξ) -b = 0, (8)
where [START_REF] Ortega | The Newton-Kantorovich theorem[END_REF] has been obtained by a change of variable u = g(ξ). Note that for technical reasons we will also extend g to the whole R N . This will be done in a way that ensures that g is convex and continuously differentiable on R N . We then show that F (u), = l, r remains concave, that F (u) exists and has a nonnegative inverse for all u ∈ R N . Therefore Newton's iterates corresponding to [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] and ( 8) converge monotonically. The numerical experiment presented in Section 3 shows that the performance of the preconditioned methods is superior compare to the original formulation of (1), or alternatively to the change-of-variable approaches.

The remainder of the article is organized as follows. In Section 2 we prove the existence and uniqueness of the solution to (1), we present the monotone convergence result for Newton's method applied to the systems (1), ( 7) and ( 8); in particular we show that Newton's method applied to [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] and ( 8) converges independently of the initial guess. In addition in Section 2.1 we deal with the fact that in practice the function g is not evaluated exactly, and we show that a two-level nested Newton's method applied to [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] still exhibits the global convergence.

Convergence analysis

In this section we analyze the convergence of Newton's method applied to the problems (1), [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] and [START_REF] Ortega | The Newton-Kantorovich theorem[END_REF]. To begin with, we present a version of the Monotone Newton Theorem and establish the existence and uniqueness of the solution of (1). Although those two results are quite standard, the proof will be presented for the reader's convenience. Then, the monotone convergence of Newton's method is established for systems (1), ( 7) and [START_REF] Ortega | The Newton-Kantorovich theorem[END_REF]. Finally in the subsection 2.1 we investigate the convergence of the preconditioned methods when g is calculated only approximatively.

The analysis presented in this section uses the notions of concavity and inverse isotonicity, so let us recall those definitions. For a more detailed discussion we refer to [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF]. Let D be an open convex subset of R N and let F : D → R N be Gâteaux differentiable. We say that F is concave if

F (u) -F (v) ≤ F (v)(u -v) (9)
for any u, v ∈ D, and we say that F is inverse isotone if

F (u) ≥ F (v) ⇒ u ≥ v (10)
for any u, v ∈ D; in addition an inverse isotone mapping F is strictly inverse isotone if [START_REF] Vázquez | The Porous Medium Equation -Mathematical theory[END_REF] holds with strict inequalities. Let us remark that inverse isotonicity implies that the equation F (u) = 0 has at most one solution. We state below a simple sufficient for the strict inverse isotonicity. For further in-depth discussion of this topic we refer to [START_REF] Rheinboldt | On M-functions and their application to nonlinear Gauss-Seidel iterations and to network flows[END_REF].

(u, v) such that J(u, v) -1 ≥ 0 and F (u) -F (v) ≤ J(u, v)(u -v). ( 11
)
Then F is strictly inverse isotone.

Let F be a Gâteaux differentiable concave mapping, such that F (u) -1 exists and is nonnegative, then, in view of (9), F satisfies the assumptions of Proposition 2.1 with J(u, v) = F (u). Similar result holds for a convex mapping with J(u, v) = F (v). On the other hand, thanks to the mean value theorem, Proposition 2.1 holds for the nonlinear mappings in the left-hand-side of (1), [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] or [START_REF] Ortega | The Newton-Kantorovich theorem[END_REF] without convexity/concavity assumption. Now, let us present a simplified version of the Monotone Newton Theorem from [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF] (theorem 13.3.4). Note that in contrast with [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF], the monotone convergence result presented below deals with concave mappings. The proof of Theorem 2.1 below is almost identical to the proof given in [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF]. Theorem 2.1 (Monotone Newton Theorem) Let D be an open convex subset of R N and let F : D → R N be continuous, Gâteaux differentiable and concave, and suppose that F (u) has a nonnegative inverse for all u ∈ D. Assume in addition that there exists u s ∈ D satisfying F (u s ) = 0 and u 0 ∈ D such that F (u 0 ) ≤ 0. Then the sequence

u n+1 = u n -F (u n ) -1 F (u n ), n ≥ 0 ( 12
)
is well defined, satisfies u n ≤ u n+1 ≤ u s and F (u n ) ≤ 0 for all n ≥ 0, and is convergent. If in addition there exists an invertible P ∈ M(N ) such that F (u n ) -1 ≥ P ≥ 0 for all n ≥ 0, then the sequence u n converges to u s .

Proof: Assume that F (u n ) ≤ 0 for some n ≥ 0 (e.g. n = 0), this implies, in view of Proposition 2.1, that u n ≤ u s . Since F (u n ) -1 is nonnegative we deduce from ( 12) that u n+1 ≥ u n . On the other hand we have that

F (u n ) -F (u s ) ≥ F (u n )(u n -u s )
which implied in view of ( 12) that

u n+1 = u n -F (u n ) -1 (F (u n ) -F (u s )) ≤ u s .
This shows that u n+1 satisfy u n ≤ u n+1 ≤ u s , and in particular F (u n+1 ) is well defined. Using concavity of F and ( 12) we have that

F (u n+1 ) -F (u n ) ≤ F (u n )(u n+1 -u n ) = -F (u n ), (13)
and, thus F (u n+1 ) ≤ 0. We have shown that the sequence (u n ) n remains in D, is nondecreasing and bounded from above, and hence converges to some u ∈ D. Let us prove that u = u s , since

F (u) -1 ≥ P we deduce that u n+1 -u n ≥ -P F (u n ) ≥ 0, implying that lim n→∞ F (u n ) = 0 since P is nonsingular.
From continuity and inverse isotonicity of F we deduce that u = u s .

Remark 2.1 Assume that F is such that u -F (u) -1 F (u) ∈ D for all u ∈ D, this is true for example if D = R N .
Then, the algorithm (12) is convergent for any initial guess; in particular the sequence (u n ) n is monotone starting from n = 1. To see that, we remark that the estimate F (u n+1 ) ≤ 0 in the proof of Theorem 2.1 resulting from (13) does not depend on the sign of F (u n ).

In fact (13) is valid as soon as u n+1 is in D.

Recall that the mappings F l and F r introduced in ( 7) and ( 8) rely on the function g which has not yet been compliantly defined. The function g coincides on R N ≥0 with the inverse of f , let us define it on the whole R N . Because the functions u → f i (u) -1 are continuous, increasing and bounded in the vicinity of zero they can be extended by continuity to u = 0. We then define g : R N → R N as a diagonal mapping, whose components g i are given by

g i (u) = f -1 i (u), u ≥ 0, f i (0) -1 u, u < 0. ( 14
)
Since f i (0) = 0, the functions g i are continuous on R; moreover, they are continuously differentiable, and, since f i are decreasing, we deduce that the functions g i are convex. Clearly the mappings F l and F r defined by ( 7) and ( 8) are concave. Let us show that the inverse of F (u), = l, r is nonnegative for all u ∈ R N . In fact the following proposition holds.

Lemma 2.1 Assume that (A 1 ) is satisfied, then the matrix F (u), = l, r is an M-matrix satisfying F (u) ≤ I ≤ F (u) -1 for all u ∈ R N ; moreover F , = l, r is strictly inverse isotone.

Proof: Let us first remark that in view of ( 14) we only need to prove the statement for u ∈ R N ≥0 . Let us denote

F u (u) = f (u) + Au -b, (15)
and w = b -Au ≥ 0. Let ε ∈ R N >0 , since g is increasing and A ≤ 0 we have that

F l (u) = I + g (w)A ≥ I + g (w + ε)A = f (g(w + ε)) -1 F u (g(w + ε)) (16)
and

F r (u) = I + Ag (u) ≥ I + Ag (u + ε) = F u (g(u + ε))f (g(u + ε)) -1 . (17)
We remark that the use of ε in the inequalities above is motivated by the fact that F u is only defined on R N >0 and not on R N ≥0 . Since f (g(•)) is a positive and diagonal matrix, we deduce from (A 1 ) that the right-hand-sides of ( 16) and ( 17) have a positive inverse. In addition, since A ≤ 0 and g is nonnegative, we deduce from Lemma 2.2 below that F (u) -1 ≥ 0 for = l, r; moreover F (u) is an M-matrix and satisfies F (u) ≤ I, which implies in turn that F (u) -1 ≥ I.

Finally, because g is diagonal and has continuously differentiable components one deduce from the mean value theorem that

F (u) -F (v) = F (z)(u -v), = l, r,
with some z ∈ R N ≥0 . In view of Proposition 2.1 this implies that F is strictly inverse isotone.

Proposition 2.2 (Existence and uniqueness of the solution) Assume that (A 1 ) is satisfied, then the solution to (1) exists and is unique.

Proof: Let us consider the mappings G l : u → g(b -Au), and let us show that G l is a contraction. Since G l ≥ 0, we have that F l (u) = I -G l (u) is a weak regular splitting of F l (u) for all u. In view of Lemma 2.1, the matrix F l (u), has a nonnegative inverse and we deduce from 2.4.17 of [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF] that ρ(G l (u)) < 1 for all u. This shows that G l is contractive on R N (with respect to some appropriate norm), and, thus G l has a unique fixed point u s ; moreover the sequence generated by

u n+1 = G l (u n ) ( 18
)
converges u s for any u 0 . Since F l (0) ≤ 0, it follows from Lemma 2.1 that u s ≥ 0, and because the restriction of g on R N ≥0 is a bijection we deduce that u s is the unique solution of (1). Let us remark that the proof of Proposition 2.2 does not rely on the concavity of F . In addition, since [START_REF] Brugnano | Iterative solution of piecewise linear systems for the numerical solution of obstacle problems[END_REF] can be expressed as

u n+1 = u n -F l (u n ) ( 19
)
we observe that the iterative Jacobi process (18) converges component-wise monotonically. We also note that, since F (u) ≤ I, we can interpret [START_REF] Casulli | A Nested Newton-Type Algorithm for Finite Volume Methods Solving Richards' Equation in Mixed Form[END_REF] as a modified Newton method, where F l (u n ) -1 has been replaced by I, which is a nonnegative subinverse of F l (u). We refer to [START_REF] Ortega | Monotone iteration for nonlinear equations with applications to Gauss-Seidel methods[END_REF] for the analysis of other Newton-like methods of this kind. Let us also note that the stationary iterations

ξ n+1 = b -Ag(ξ n )
corresponding to the system (8) converge to ξ s = f (u s ).

We now in the position to prove that Newton's method applied both to the original problem formulation (1) and the preconditioned problems ( 7) and ( 8) converges monotonically. Remark however that, since f may be unbounded at the origin, the mapping F u from ( 15) is only well defined on D = R N >0 . On the other hand the initial guess u 0 required by Theorem 2.1 needs to satisfy u 0 ≤ u s , while u s does not have to be strictly positive. Therefore, unless some additional hypotheses are made, Newton's method may be inapplicable to the original formulation (1). Proposition 2.3 (Convergence of the original method) Assume that b > 0, then there exists an initial guess u 0 > 0 such that Newton's method applied to (1) converges monotonically.

Proof: Let F u be given by [START_REF] Pao | Accelerated monotone iterative methods for finite difference equations of reactiondiffusion[END_REF] and let D = R N >0 , in view of the assumption (A 1 ) the mapping F u defined on D is continuously differentiable and concave, in addition F u (u) has a nonnegative inverse for all u ∈ D. It remains to show that u s ∈ D and that there exists u 0 ∈ D satisfying F u (u 0 ) ≤ 0.

Let 1 N denote the element of R N with all unit components, from continuity of f and the fact that f (0) = 0 we deduce that there exists > 0 such that f (1 N ) ≤ b, and therefore F u ( 1 N ) ≤ 0. This implies that u 0 = 1 N > 0 is an appropriate initial guess. Proposition 2.4 (Convergence of the preconditioned methods) Newton's method applied to [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] and (8) converges for any initial guess. In particular the sequence of Newton iterates (u n ) n converges monotonically starting from n = 1.

Proof: It follows form Lemma 2.1 and Proposition 2.2 that F , = l, r satisfies the assumptions of Theorem 2.1 with D = R N and u 0 = 0. The global convergence follows from Remark 2.1.

Remark 2.2

In order to fit the problems [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] and (8) into the framework of Theorem 2.1 we have defined the mapping g is an extension of f -1 to the whole R N . This is however a rather theoretical construction, since in view of Proposition 2.4 the iterates starting from any u 0 ≥ 0 such that F (u 0 ) ≤ 0, = l, r (e.g. u 0 = 0) will remain in R N ≥0 .

Remark 2.3 Let us note that the convergence analysis presented above applies to some mildly nonlinear systems that can not be written in the form of (1). Let us consider the system

F (τ ) := v(τ ) + Lu(τ ) -b = 0 ( 20
)
where L ∈ M(N ), while v and u are the diagonal mappings from R N to R N that are nondecreasing, but not necessarily strictly increasing. The system (20) typically results from the discretization of some constraint PDEs. Examples of problems leading to (20) include degenerate Richards' equation [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF], the evolutionary dam problem [START_REF] Menéndez | On the uniqueness of the solution of the evolution dam problem[END_REF], Stefan or Hele-Shaw problems [START_REF] Vázquez | The Porous Medium Equation -Mathematical theory[END_REF], as well as some classical elliptic or parabolic obstacle problems [START_REF] Brugnano | Iterative solution of piecewise linear systems for the numerical solution of obstacle problems[END_REF]. Let D be the diagonal of L and A = D -L, then, denoting ψ(u) = v(τ ) + Du(τ ), we can express the system (20) as

ψ(τ ) + Au(τ ) -b = 0.
Assume that ψ is strictly increasing, then, using a new variable ξ = ψ(τ ), and denoting g(ξ) = u(ψ -1 (ξ)), we obtain the system

ξ + Ag(ξ) = b (21)
similar to [START_REF] Ortega | The Newton-Kantorovich theorem[END_REF]. Now, if u is merely nondecreasing, then g is not bijective and (21) can not be cast into (1). Nevertheless preconditioned system similar to (7) can be obtained in the following form:

Find ξ such that ξ = b -Au with u -g(b -Au) = 0. (22)
It is easy to show that ξ is a solution of [START_REF] Dolean | Nonlinear preconditioning: how to use a nonlinear Schwarz method to precondition Newton's method[END_REF] if and only if it solves [START_REF] Menéndez | On the uniqueness of the solution of the evolution dam problem[END_REF]. Now, let us show that the systems (21) and ( 22) can be fitted into the framework of the Monotone Newton Theorem. Assume that F form (20) is defined on R N and that F (τ ) is an M-matrix for all τ , then on shows that the I -g(b -Au) and I + Ag(ξ) are also M-matrices whose inverses are bounded from below by I. Assume in addition that F (τ ) = 0 has a solution, then in order to apply Theorem 2.1 it remains to show that g i are concave for all i ∈ {1, . . . N }. In order to do so let us assume that v is concave and u is convex. Denoting ζ = ψ -1 i (ξ), we have

g i (ξ) = u i (ζ)ψ i (ζ) -1 = u i (ζ) v i (ζ) + D i u i (ζ)
.

Since the function γ(p, q) = p q + D i p , q, p ≥ 0 is nonincreasing with respect to q and nondecreasing with respect to p we deduce that

g = ∂γ ∂p u i + ∂γ ∂q v i (ψ -1 i ) is nonnegative.
The numerical experiment presented in Section 3 provides the evidences that the preconditioning substantially improves the convergence of Newton's method. To support this observation theoretically we present the following proposition stating that the preconditioned methods lead to a larger solution updates. Proposition 2.5 Let u ∈ R N ≥0 be such that F u (u) ≤ 0 and f (u) < +∞. Let u up , u l up and u r up denote the update generated by Newton's method applied to (1), ( 7) and (8) respectively, starting from the initial guess u. Then, u l up ≥ u up and u r up ≥ u up . and there exists a sequence (σ n ) n ≥ 0 such that lim n→∞ σ n = 0 and

-σ n ≤ F (u n ) -F n, , (28)
then u n converges to u s .

Proof: Let F (u n ) ≤ 0 for some n ≥ 0 (e.g. for n = 0), since F n, ≤ 0 we deduce from (24) that u n+1 ≥ u n . Let us show that u n+1 ≤ u s . From ( 24), ( 25), (26) and using concavity of F we deduce that

u n+1 ≤ u n -F (u n ) -1 F (u n ) = u n -F (u n ) -1 (F (u n ) -F (u n )) ≤ u s .
This implies in particular that u n+1 ∈ D. It follows from concavity of F that

F (u n+1 ) -F (u n ) ≤ F (u n )(u n+1 -u n ),
and using (24) we obtain

F (u n+1 ) -F (u n ) ≤ -F (u n )J -1 n, F n, or F (u n+1 ) ≤ I -F (u n )J -1 n, F n, + F (u n ) -F n, .
Since J -1 n, ≥ 0 we deduce from (26) that

I -F (u n )J -1 n, ≥ 0,
and, in view of (25), we deduce that F (u n+1 ) ≤ 0. The sequence (u n ) n is nondecreasing and bounded from above; therefore (u n ) n converges to some u. Now, assume that ( 27) and (28) are satisfied. Combining ( 27) and ( 24) we find that 0 ≥ -P lim n→∞ F n, ≥ 0, which implies that lim n→∞ F n, = 0. In turn, the condition (28) and the continuity of F yield F ( u) = 0, and in view of Proposition 2.1, we deduce that u = u s .

To complete this section we show that the nested Newton's method applied to the problem (7) satisfies the assumptions of Proposition 2.6. We begin with a following technical lemma.

B 1 = D 1 -M 1 ≥ 0, B 2 = D 2 -M 2 ≥ 0. Since, M 2 ≥ M 1 we deduce that D 2 ≥ D 1 ≥ 0 and B 1 ≥ B 2 ≥ 0, and therefore 0 ≤ D -1 2 B 2 ≤ D -1 1 B 1 .
Denoting by ρ(M ) denote the spectral radius of a matrix M , we deduce from 2.4.8 and 2.4.17 of [START_REF] Ortega | Iterative Solutions of Nonlinear Equations in Several Variables[END_REF] that ρ(D -1 2 B 2 ) ≤ (D -1 1 B 1 ) < 1, and that M -1 2 ≥ 0. The following proposition draws the connection between Proposition 2.6 and the approximate evaluation of the function g.

= I + f (w ) -1 A, then F l (u) ≤ F ≤ 0, (29)
and J ≥ F l (u) and J -1 ≥ I.

β(u n i ) + ∆t h 2 j∈Ni (u n i -u n j ) = β(u n-1 i ) + ∆tqδ i,1 , (34)
where δ i,1 stands for the Kronecker symbol and where N i denotes the set of cells neighboring with

K i N i =    {2} i = 1, {i -1, i + 1} 1 < i < N, {N -1} i = N.
Let L denote the tridiagonal matrix associated to the discretization of the diffusion operator in the left-hand-side of (34) and let, for a given n ≥ 0, b n denote the right-hand-side of (34). Imposing (34) for all i ∈ {1, . . . , N } leads to the following system of algebraic equations

β(u) + ∆t h 2 u + L - ∆t h 2 I u = b n , (35)
which has to be solved for each n ∈ {1, . . . , N T }. It is easy to verify that b n ≥ 0 and that

f (u) = β(u) + ∆t h 2 u and A = L - ∆t h 2 I satisfy the assumptions (A 1 ) -(A 2
). The objective of the numerical experiment is to evaluate the efficiency of Newton's method (NM) applied to left and right-preconditioned problems

F n l (u) := u -g(b n -Au) = 0 (36) and F n r (u) := u + Ag(u) -b n = 0. ( 37
)
Those Jacobi-Newton methods are compared, in terms of the performance, with three more traditional approaches specified below. u-formulation: NM applied to (35) in the original form

F n u (u) := β(u) + Lu -b n = 0. ( 38
)
In view of Proposition 2.3 this method is monotonically convergent provided that the initial guess satisfies F (u 0 ) ≤ 0. and

err u, m, = u •, m, -u •, m,ref L ∞ (0,T ;L 1 (0,1)) u •, m,ref L ∞ (0,T ;L 1 (0,1))
Because the qualitative behaviour of err u, m, of err v, m, is similar we will be only reporting the latter error metric. Performance comparison. The first set of tests is performed using the fixed mesh size parameter N = 100. In accordance with the results reported in [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF], Figures 2(a) and 2(b) witness the qualitative differences in the performance of u and v-formulations on one hand, with v-formulation being several times faster, and the performance of v and τ -formulations on the other hand, with τ -formulation providing the most significant speedup. It can also be noted that, in contrast with u and v-formulations, the sensibility of the τ -formulation to the parameter m is very limited. In turn, Figure 2(c) shows a relatively similar behaviour of τ -formulation and the preconditioned methods, with the latter ones requiring an even smaller number of iterations. Computational overhead due to local problem solution. As it can be observed on Figure 2(c) the preconditioned Newton's methods require fewer iterations than the method based on τ -formulation. However, each iteration of the Jacobi-Newton method requires to solving a set of the scalar nonlinear equations. Those inner calculations produce a certain computational overhead. To access the overall computational effort required by the preconditioned methods we present the analysis in terms of the CPU time. Figures 3(a

Conclusion

For the problems involving only diagonal nonlinearities and satisfying Monotone Newton Theorem, we have proposed a nonlinear preconditioning procedure based on the Jacobi method. This preconditioning is computationally inexpensive and leads to a monotone Newton's method that converges globally and faster than the original one. We believe that this method is particularly efficient for problems involving stiff nonlinearities. This point is illustrated by the numerical experiment based on the porous media equation. We observe that the convergence of the original Newton's method is very slow and deteriorates as the diagonal nonlinearity gets stiffer. In contrast, our newly proposed method exhibits a fast convergence independently of the nonlinear stiffness, which in some sense is absorbed by the preconditioner. The preconditioned method also turns out to be more efficient than the alternative nonmonotone methods based on the change of the variable.

Proposition 2 . 1 (

 21 Inverse isotonicity) Let D be an open convex subset of R N and let F : D → R N , suppose that for any u, v ∈ D there exists a nonsingular matrix J

Lemma 2 . 2 2 ≥ 0 .

 2220 Let M 1 , M 2 ∈ M(N ) with M 1 being an M-matrix, and M 2 having nonpositive offdiagonal elements and satisfying M 2 ≥ M 1 , then M -1 Proof: Let D 1 and D 2 denote the diagonal of M 1 and M 2 respectively and

Proposition 2 . 7

 27 Let u ∈ R N ≥0 be such that F l (u) ≤ 0. Let w denote the unique solution of f (w) = b -Au and let w satisfy u ≤ w ≤ w. Let F = u -w and J

Figure 1 :

 1 Figure 1: Approximate solution at different time steps

  (a) err v, m, for v-formulation (solid blue) and u-formulation (dashed blue) (b) err v, m, for v-formulation (solid blue) and τ -formulation (magenta) (c) err v, m, for τ -formulation (magenta), right preconditioned (black) and left preconditioned (red) Newton's method

Figure 2 :

 2 Figure 2: Relative error err v, m, as the function of the average number of Newton's iterations per time step

  ) and 3(b) show, for different values of the mesh size parameter N ∈ {200, 400, 800, 1200}, the comparison of left (respectively right) preconditioned NM with the method based on τ -formulation. It can be observed that for all except very small prob-

  (a) err v, m, for left preconditioned NM (solid lines) and τ -formulation (dashed lines) (b) err v, m, for right preconditioned NM (solid lines) and τ -formulation (dashed lines)

Figure 3 :

 3 Figure 3: Relative error err v, m, as the function of CPU time for different grid sizes

Proof: Let us first consider the system [START_REF] Kantorovich | On Newtons method for functional equations[END_REF], and let us denote w = b -Au. Since F u (u) ≤ 0 we deduce that f (u) ≤ w and, thanks to the mean value theorem, we have that

for some z satisfying f (u) ≤ z ≤ w. One the other hand, F l (u) = I + g (w)A = I + f (g(w)) -1 A.

Therefore u l up satisfies the equation f (g(z)) I + f (g(w)) -1 A u l up -u = -F u (u), while u up satisfies f (u)

Since f is nonincreasing, A is nonpositive and u ≤ g(z) ≤ g(w) we have that

In view of (A 1 ) both sides of the above inequality are the M-matrices, therefore, we deduce that u l up ≥ u. Now, we consider the system (8) and we denote ξ = f (u) and ξ up = f (u r up ). Writing down a single step of Newton's method, and using again the mean value theorem, we have

for some z satisfying u ≤ z ≤ u up . In view of [START_REF] Cai | Nonlinearly preconditioned inexact Newton algorithms[END_REF] and observing that f (z) ≤ f (u) and deduce that u r up ≥ u up .

Convergence of the inexact methods

The application of Newton's method to the preconditioned problems [START_REF] Kantorovich | On Newtons method for functional equations[END_REF] and ( 8) requires evaluation of the function g, which in general can not be done exactly. In order to compute g(v) for some v ∈ R N ≥0 one has to solve a set of scalar nonlinear equations of the form f (w) = v. This can be achieved by any appropriate iterative method, such as bisection, regula falsi or Newton's method again. The fact that in practice the function g is evaluated only approximatively gives rise to the following sequence of the inexact iterations

Here F n, and J n, denote some approximations of F (u n ) and F (u n ) respectively. Let us give the conditions under which the inexact method (24) converges to u s .

Proposition 2.6 Let D be an open convex subset of R N and let F : D → R N be continuous, Gâteaux differentiable and concave, and suppose that F (u) have a nonnegative inverse for all u ∈ D. Assume in addition that there exists u s ∈ D satisfying F (u s ) = 0 and u 0 ∈ D such that F (u 0 ) ≤ 0.

Let (u n ) n be a sequence constructed by the following algorithm: For all n ≥ 0

3. Use (24) to compute u n+1 .

Then, the sequence (u n ) n is well defined for all n ≥ 1 and satisfy u n ≤ u n+1 ≤ u s and F (u n ) ≤ 0 for all n ≥ 0.

If in addition

there exists an invertible P ∈ M(N ) such that J -1 n, ≥ P ≥ 0 for all n (27)

Proof: Note that

and

Subtracting (32) from (31) we find

which, combined with (32), implies (29). Since w ≤ w and since f is diagonal and concave we have that

which implies that J ≥ F l (u). In turn it follows from Lemma 2.2 that J -1 ≥ 0, and therefore

Let us show that, for a given u n satisfying F l (u n ) ≤ 0, the computation of an approximation w n, of w n satisfying Proposition 2.7 can be achieved by Newton's method. Let r n = b -Au n , since F l (u n ) ≤ 0 and in view of (31) we have that

Using similar arguments as in the proof of Theorem 2.1 one shows that u n ≤ w n,k ≤ w n,k+1 ≤ w n for all k ≥ 0 and that the sequence (w n,k ) k converges toward w n . In view of Proposition 2.7 we have that for any k ≥ 0 the quantities

satisfy (25), ( 26) and ( 27). The extraction of the sequence w n, providing (28) can be done by setting w n, = w n,κ(n) where κ(n) is the smallest integer satisfying w n -w n,κ(n) ≤ 10 -n .

Remark 2.4

The result similar to Proposition 2.7 can be established for the right-preconditioned method. In that case one has to require that w satisfies w ≤ w and u ≤ b -Aw .

However, in contrast with the left preconditioned method, it is unclear how such approximated values w can be constructed in practice.

Numerical experiment

Let us consider the porous medium equation (see [START_REF] Vázquez | The Porous Medium Equation -Mathematical theory[END_REF])

on (0, 1) × (0, T ). The nonlinearity in the accumulation term is given by β(u) = u 1/m with m > 1.

We consider the Neumann boundary conditions

with q > 0, and the constant initial condition u(x, 0) = u ini > 0. The value of u 0 is going to be chosen close to zero leading to a traveling wave like solution. Figure 1 exhibit the approximate profile of β(u) at different time steps computed for m = 10, q = 10 4 and the time step ∆t = 1.2 10 -4 . Equation (33) is discretized using the standard finite volume method with the time integration performed by the backward Euler scheme.

Let N be a positive integer, let h = 1/N and let x i = i/N for i ∈ {0, . . . , N }; the set of the finite volumes (K i ) i∈{1,...,N } is defined by K i = (x i-1 , x i ). Let N T be a positive integer and let v-formulation: The problem (35) is reformulated with respect to the variable v with u = β -1 (v) and NM is applied to

τ -formulation: Following [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF] we introduce the function pair τ → (u(τ ), v(τ )) such that

for all τ and max (u (τ ), v (τ )) = 1.

Then NM is applied to

At each time step n and for each of the formulations (36)-( 40) the sequence of the approximate solutions is computed using Newton's method

until the stopping criterion

is satisfied for some small predefined tolerance parameter . As the initial guess we use the value of the variable obtained at the previous time step (this value will obviously differ from one formulation to another). This choice of the initial guess is motivated by the following observation.

Remark 3.1 The solution of (33) (under given initial and boundary conditions) satisfies ∂ t u ≥ 0. This property is reproduced at the discrete level by the approximate solution resulting from u-formulation and the preconditioned methods. For = u, r, l, let us denote by u n the approximate solution of F n (u) = 0, then one can show that F n (u n-1 ) ≤ 0, and therefore u n 0 = u n-1 provides an appropriate choice of the initial guess. Let us give the proof by induction for the case of u-formulation. The proof for preconditioned methods is similar, given that for all u ∈ R N ≥0 , all n and = r, l one has

u (u 0 ) < 0 providing, in view of Proposition 2.3, that the sequence of Newton's iterates are monotonically increasing, and that u 1 satisfies u 1 ≥ u 0 and F 1 u (u 1 ) ≤ 0. Next, we show that if the statement

is true for some n = p ≥ 1, then it is true for n = p + 1. To do that we notice that for n ≥ 1

Therefore, if (41) is satisfied for some n = p ≥ 1, then F p+1 u (u p ) ≤ 0, which implies u p+1 ≥ u p and F p+1 u (u p+1 ) ≤ 0 in view of Proposition 2.3. Now, we present the results of the numerical experiment. The test case is configurated as follows: in order to allow for the use of u-formulation we chose a positive initial condition β(u ini ) = 10 -10 , we set q = 10 4 , T = 1.2 10 -2 , N T = 100 and we let the parameter m to take values in the set {4, 8, 16, 32}. For a given value of m, the tolerance and a specific solution method we denote by u n, m, n∈{1,...,N T } and v n, m, n∈{1,...,N T } the approximate solution of (35). The methodology of the numerical experiment is similar to [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF], that is for each value of m we compute, using τ -formulation and the tolerance ref .


Nazli Mohammad

Factors Influencing the Development of Trust

Introduction

Scholars in various disciplines have considered the causes, nature, and effects of trust. In this research, prior approaches to studying trust are considered, including relation between belonging with a community and trust, and social characteristics. This research addresses cultural factors that may have a direct effect on the level of trust in society. Specific questions to be addressed are:

What is trust? What are factors influencing the development of trust?

Research Goal

Understanding relationship between individual values, big five personalities, social isolation, and anomie with trust. The purpose is to illuminate the key role of trust in critical social processes.

Methods

The data for this study were drawn from the World Values Survey, PWC's 2016 Global CEO Survey, an independent survey firm in February 2016 from a nationally representative sample of 1,095 working adults in the U.S., and an independent survey in 2019 from 127 students of the MBA at the University of Nevada, Reno.

Definition of trust

There are endless numbers of trust definitions. But my favorite one is that "I am OK with feeling vulnerable, giving you the opportunity to harm me but I believe that you will not ".

Generally speaking, would you say that most people can be trusted? Who are in your circle of trust?

Determinants of trust

What are the factors influencing the development of trust?

First individual values, based on the Schwartz value survey there are 10 personal values like selfenhancement, power, achievement, openness to change. Our study shows that these values can positively affect trust.

We also found that personality can have an impact on trust through social isolation. The big 5 personality test measures extraversion, agreeableness, openness, conscientiousness, and neuroticism. Respondents who score high on extraversion are less likely to indicate social isolation. On the other hand, those who score high on Neuroticism are more likely to indicate social isolation. To a lesser extent, agreeableness and Openness are negatively related to social isolation.

Perceived social isolation is also known as loneliness. More isolation could result in a lower level of trust. For example, respondents who score high on extraversion are less likely to indicate social isolation. On the other hand, those who score high on Neuroticism are more likely to indicate social isolation. We could not find significant correlations between PSI and Cultural/Individual Values. Anomie captures the extent to which people feel that their community has rules that are followed by all. Questions like you can break any rules if they do not seem to apply or I follow whatever rules I want to follow, show a mismatch between the group and individual standards. People feeling disconnected from society that results in less trust.

Conclusion

Trust is measured as investment, investment of time, effort, or money. It becomes evident that without trust, cooperation, interactions, and trust-based relationships could never be established. So, taking together we need trust. Trust is indispensable. Trust empowers us and we can empower others by trusting them.