Reza Mollapourasl

Ali Fereshtian
email: a.fereshtian@sru.ac.ir

Florin Avram
email: florin.avram@univ-pau.fr

Mohan Kadalbajoo

Mehdi Dehghan

RBF approximation by partition of unity for valuation of options under exponential Lévy processes

Keywords: Radial basis functions, partition of unity, Option pricing, L\'{e}vy processes, Markov processes Radial basis functions, partition of unity, Option pricing, Lévy processes, Markov processes

The prices of some European and American-style contracts on assets driven by a class of Markov processes containing, in particular, L\'{e}vy processes of pure jump type with infinite jump activity, are obtained numerically, as solutions of the partial integro-differential equations (PIDEs) they satisfy. This paper overcomes the ill-conditioning inherent in global meshfree methods by using localized RBF approximations known as the RBF partition of unity (RBF-PU) method for (PIDEs) arising in option pricing problems in L\'{e}vy driven assets. Then, Crank-Nicolson, Leap-Frog (CNLF) is applied for time discretization. We treat the local term using an implicit step, and the nonlocal term using an explicit step, to avoid the inversion of the nonsparse matrix. For dealing with early exercise feature of American option and solving free boundary problem we use the implicit-explicit method combined with a penalty method. Efficiency and practical performance are demonstrated by numerical experiments for pricing European and American contracts.

Dear Editor,

We would like to submit our manuscript entitled

RBF approximation by partition of unity for valuation of options under exponential Lévy processes

for exclusive consideration of publication as a research paper in Journal of Computational Science.

The prices of some European and American-style contracts on assets driven by a class of Markov processes containing, in particular, Lévy processes of pure jump type with infinite jump activity, are obtained numerically, as solutions of the partial integro-differential equations (PIDEs) they satisfy. This paper overcomes the ill-conditioning inherent in global meshfree methods by using localized RBF approximations known as the RBF partition of unity (RBF-PU) method for (PIDEs) arising in option pricing problems in Lévy driven assets. Then, Crank-Nicolson, Leap-Frog (CNLF) is applied for time discretization. We treat the local term using an implicit step, and the nonlocal term using an explicit step, to avoid the inversion of the nonsparse matrix. For dealing with early exercise feature of American option and solving free boundary problem we use the implicit-explicit method combined with a penalty method. Efficiency and practical performance are demonstrated by numerical experiments for pricing European and American contracts.

Thank you for considering our paper. We appreciate your effort and look forward to your response.

With best regards

Ali Fereshtian Reza Mollapourasl (corresponding author, mollapor@oregonstate.edu) Florin Avram

We proposed RBF-PU method for spatial discretization of PIDE to price American and European options under Lévy model. Also, penalty method is used to solve the free boundary problem arisen in pricing American option. Then, CNLF is applied for time discretisations and we treat the local term using an implicit step and the nonlocal term using an explicit step to avoid to avoid the inversion of the nonsparse matrix. These result in a linear algebraic system with a sparse matrix that has a small condition number. The shape parameter in the RBF affects the accuracy and stability of the numerical methods. Numerical results confirm that RBF-PU method is less sensitive to the change of shape parameter. An increase of the number of nodal points and correspondingly an increase of the number of patches in RBF-PU method also leads to an improvement of the approximation. The effect of the time discretization is measured by studying the temporal error. For the American and European option cases where the CNLF scheme is combined with penalty method, we conclude that the rate of convergence is of at least first order in time for more cases. Also, numerical results confirm that our scheme is stable and second order convergent in space. The experiments also demonstrated that European and American option prices with error around 1.0e -4 can be computed in less than one second on a PC. Thus, the developed method is very fast and accurate.

Introduction

It is widely recognized that the classic option pricing model proposed in 1973 by Black and Scholes in [6] and Merton in [START_REF] Merton | Option pricing when underlying stock returns are discontinuous[END_REF], does not ideally fit observed empirical market data. Two identified empirical features have been under much attention, the first one is skewed distribution with higher peak and heavier tails of the return distribution and the second one is the volatility smile [START_REF] Andersen | Jump-diffusion processes: Volatility smile fitting and numerical methods for option pricing[END_REF]. Jumps are regularly observed in the discrete movement of stock price and these jumps can not be captured by the log normal distribution characteristic of the stock price in the Black Scholes model. Therefore an alternative model is necessary to overcome these issued. To resolve these, issue several models have been proposed in the literature. Among these, the jump diffusion model introduced by Merton [35] and Kou [START_REF] Kou | A jump-diffusion model for option pricing[END_REF] is one of the most used model. Merton proposed a log-normally distributed process for the jump-amplitudes, while Kou suggested logdouble-exponentially distributed process.

These models have finite jump activity, unlike the more general approach with possibly infinite jump activity -see [START_REF] Carr | The fine structure of asset returns: An empirical investigation[END_REF][START_REF] Cont | Financial Modelling with Jump Processes[END_REF]. Another motivation for our study comes from stochastic volatility models with jumps. The model proposed by Bates in [START_REF] Bates | Jumps and stochastic volatility: Exchange rate processes implicit in Deutsche mark options[END_REF] with jumps only in the value of the underlying asset is an example of such an approach. More general jump-diffusion models with stochastic volatility are considered in [START_REF] Eraker | The impact of jumps in volatility and returns[END_REF], for example.

A variety of numerical techniques have been proposed for pricing options under both models with stochastic volatility and models with jumps. In particular, in the case of models with jumps, we mention the paper by Matache et al. [START_REF] Matache | Wavelet galerkin pricing of American options on Lévy driven assets[END_REF], where a finite element method is developed, the papers by Andersen and Andreasen [START_REF] Andersen | Jump-diffusion processes: Volatility smile fitting and numerical methods for option pricing[END_REF], by Cont and Voltchkova [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF], by d'Halluin et al. [START_REF] Halluin | A penalty method for American options with jump diffusion processes[END_REF][START_REF] Zvan | Penalty methods for American options with stochastic volatility[END_REF], and by Clift and Forsyth [START_REF] Clift | Numerical solution two asset jump diffusion models option valuation[END_REF], where finite difference schemes are presented, and the papers by Fang and Oosterlee [START_REF] Fang | A novel pricing method for European options based on fouriercosine series expansions[END_REF][START_REF] Fang | Pricing early-exercise and discrete barrier options by fouriercosine series expansions[END_REF], where a numerical method based on the characteristic function is proposed.

In models with jumps, the non local integral term leads to dense matrices after discretization [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF] and efficient numerical methods are required for pricing of complex contracts and for calibration of model parameters. This fact means that fully implicit schemes encounter difficulties such as dense matrix inversions whereas fully explicit schemes impose stability restrictions [START_REF] Briani | Implicit-explicit numerical schemes for jump-diffusion processes[END_REF]. For some solutions to this problem in [START_REF] Halluin | Robust numerical methods for contingent claims under jump diffusion processes[END_REF][START_REF] Lord | A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes[END_REF] authors have used an iterative procedure for solving the discretized equations, and more efficient approaches using implicit-explicit Runge-Kutta schemes which treat the integral term explicitly were later proposed by Cont and Voltchkova in [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF] and by Briani et al. [START_REF] Briani | Implicit-explicit numerical schemes for jump-diffusion processes[END_REF]. Also, an unconditionally stable alternating direction implicit (ADI) method was proposed by Andersen and Andreasen in [START_REF] Andersen | Jump-diffusion processes: Volatility smile fitting and numerical methods for option pricing[END_REF].

Option pricing using RBFs was explored in one dimension for European and American options by Hon et al. [START_REF] Hon | A radial basis function method for solving options pricing model[END_REF], [START_REF] Wu | Convergence error estimate in solving free boundary diffusion problem by radial basis functions method[END_REF]. Pettersson et al. [START_REF] Pettersson | Improved radial basis function methods for multi-dimensional option pricing[END_REF] presented an RBF based method for multi-dimensional European options, and American options in both one and two dimensions are investigated by Fasshauer et al. [26]. All of these papers employ global RBF collocation methods. Global RBF collocation methods lead to dense linear systems, and computational costs that become prohibitive as the number of dimensions increase. In order to address the computational cost issues of the global RBF method, we need to introduce locality. An easy way to do that is to use compactly supported RBFs, such as the Wendland functions [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions minimal degree[END_REF], but then the spectral convergence properties are lost. Here we take another approach, where the infinitely smooth RBFs are still used in the approximation but over local subregions of the computational domain [START_REF] Mollapourasl | Radial basis functions with partition of unity method for American options with stochastic Computational Economics[END_REF]. The possibility of using RBFs in a partition of unity scheme was mentioned in [START_REF] Babuška | The partition of unity method[END_REF], further discussed in [START_REF] Wendland | Fast evaluation of radial basis functions: Methods based on partition of unity[END_REF], implemented for interpolation on the sphere, in the plane, and in three-dimensional space in [START_REF] Cavoretto | Spherical interpolation using the partition of unity method: efficient and flexible algorithm[END_REF][START_REF] Cavoretto | A meshless interpolation algorithm using a cell-based searching procedure[END_REF]. Also, we applied localized partition of unity method for pricing American and European option prices as well as Greeks in [START_REF] Mollapourasl | RBF-PU method for pricing options under the jump-diffusion model with local[END_REF].

We introduce below a local method for option pricing problems under Lévy processes, based on radial basis function partition of unity (RBF-PU) approximation. Some advantages of partition of unity method based on radial basis function are the flexibility in choosing ansatz spaces and their good approximation properties. Also, this method helps us to have an sparse system with low condition number and control the computational cost and run time of the algorithm. The outline of the paper is as follows.

In the next section, we describe Lévy's model and partial integro-differential equations (PIDEs) for pricing European options. Also, American option is formulated as linear complimentary problem (LCP) then we apply penalty method to solve it. In section 3, the basics of radial basis functions as well as the partition of unity method based on RBFs and the definition of the differentiation matrices are introduced. In section 4, RBF-PU method is used to discetize differential operator and for integral operator, we use trapezoid rule. Then, (CNLF) is applied for time discretization and we treat the local term using an implicit step and the nonlocal term using an explicit step to avoid the inversion of the nonsparse matrix. Finally, in the last section, the accuracy and efficiency of the proposed method is numerically investigated for European and American options.

Option Pricing under the Exponential Lévy Model

Observation of sudden, large movements in the prices of financial assets has led to the use of stochastic processes with discontinuous trajectories -jump processes -as models for financial assets. Exponential Lévy models [START_REF] Cont | Financial Modelling with Jump Processes[END_REF] offer analytically tractable examples of positive processes with jumps which are simple enough to allow a detailed study, and yet complex enough to allow calibrating the model to market prices of options. The Markov property of the price implies that option prices are solutions of partial integro-differential equations (PIDEs) which involve, in addition to a (possibly degenerate) second-order differential operator, a non-local integral term which requires specific treatment both at the theoretical and numerical level.

We consider here the class of models where the risk neutral dynamics of the underlying asset is given by S t = S 0 exp((rq)t + X t ) where X t is a Lévy process with a Lévy triplet (σ 2 , γ, ν). The law of Lévy process is completely determined by its characteristic triplet (σ 2 , γ, ν), where σ 2 is called Gaussian variance, γ is related to the drift, and the Lévy measure satisfies the following condition

R min(1, x 2 )ν(dx) < ∞.
The probabilistic interpretation of ν is that ν(dx) gives the expected number of jumps with size between x and x + dx, which the process makes between time 0 and 1. An infinite activity Lévy process is described by an infinite number of jumps, that is, ν(R) = ∞.

Example 1. The CGMY process [START_REF] Carr | The fine structure of asset returns: An empirical investigation[END_REF] is an infinite activity processes defined by the Levy measure

ν(dx) =            Ce -G|x| |x| 1+Y dx, x < 0 Ce -M x x 1+Y dx, x > 0 (1) for constants C > 0, G ≥ 0, M ≥ 0, and 0 < Y < 2.
Example 2. In the case Y = 0, a special choice of the exponents G, M yields the subclass of Variance-Gamma or VG processes, which were introduced earlier by Madan and Seneta (1990) [33] as time changed Brownian motions.

The standard tool to assign prices to options under Lévy processes consists in changing to a convenient probability measure and taking the expectation of the discounted prices and define the pricing function v driven by X t by

v(τ, x) = e -rτ E[g(x + (r -q)τ + X τ )] (2)
where g(•) is the payoff function, τ = Tt is the time to maturity T , and x = log( S K ) is the log price of the underlying asset S.

Replace now the infinitely many jumps smaller than ǫ by a Brownian motion, leaving only the jumps whose magnitudes are larger than ǫ, which are driven by a finite activity Levy model [START_REF] Asmussen | Approximations of small jumps of Lévy processes with a view towards simulation[END_REF]. To be specific, we define its corresponding Levy process X t that has the Levy triplet (σ 2 , γ, ν) given by

σ2 = σ 2 + ǫ -ǫ z 2 ν(dz) (3)
ν = ν1 |x|≥ǫ γ = - 1 2 σ 2 - R (e x -1 -x1 |x|≤1 )ν(dx)
where 1 is the indicator function. Now let

u(τ, x) = e -rτ E[g(x + (r -q)τ + Xτ )]
where u(τ, x) satisfies the following partial integro-differential equation(PIDE)

∂ τ u -Pu = 0, τ ∈ (0, T ]; u(0, x) = g(x), x ∈ R (4)
where the integro-differential operator P is

Pu(τ, x) = 1 2 σ2 ∂ 2 x u(τ, x)+(r -q - σ2 2 -λ(ǫ)ξ(ǫ))∂ x u(τ, x)-(r +λ(ǫ))u(τ, x)+ |y|≥ǫ u(τ, x+y)ν(dy)
where λ(ǫ) = |z|≥ǫ ν(dz), ξ(ǫ) = 1 λ(ǫ) |z|≥ǫ (e z -1)ν(dz).

Theorem 1. Let g be Lipschitz: |g(x)g(y)| ≤ K|x -y|. Let v(τ, x) and u(τ, x) be defined by (2) and (4), respectively. Then

|v(τ, x) -u(τ, x)| ≤ C |y|≤ǫ |y| 3 (dy) σ2 . (5)
Proof. [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF] Remark 1.

If lim x→0 ν(x)|x| 1+β = a > 0 with 0 ≤ β < 2, then (5) gives |v(τ, x) -u(τ, x)| ≤ C(β)ǫ
so the approximation error is proportional to ǫ. This case includes all practical examples used in option pricing such as variance Gamma, NIG, and tempered stable processes [START_REF] Cont | A finite difference scheme for option pricing in jump diffusion and exponential Lévy models[END_REF].

Example 3. Two important examples of payoff functions are

g(x) = max(K -Ke x , 0), for put option; max(Ke x -K, 0), for call option. (6)
In the remaining of this paper we will focus on put options, and as additional boundary conditions we impose lim

x→-∞ u(τ, x) = Ke -rτ , lim x→∞ u(τ, x) = 0. ( 7
)
For American options, it is well known that for each time τ ∈ [0, T ], there always exists a stock price x for which early exertion before final time T is advantageous [START_REF] Seydel | for Computational Finance[END_REF]. Therefore the American option pricing problem can be formulated in the form of an inequality by adopting the linear complementarity problem (LCP) formulation

   ∂ τ u(τ, x) -Pu(τ, x) ≥ 0, u(0, x) ≤ g(x), (∂ τ u(τ, x) -Pu(τ, x))(u(0, x) -g(x)) = 0, (8)
for (τ, x) ∈ R × (0, T ] and we need to impose the following boundary conditions

lim x→-∞ u(τ, x) = K, lim x→∞ u(τ, x) = 0. ( 9
)
For solving the free boundary problem, there are some techniques such as linear programming [START_REF] Borici | Pricing American put options by fast solutions of the linear complementarity problem[END_REF] or penalty methods presented in [START_REF] Zvan | Penalty methods for American options with stochastic volatility[END_REF][START_REF] Safdari-Vaighani | A radial basis function partition of unity collocation method for convection-diffusion equations arising in financial applications[END_REF][START_REF] Shcherbakov | Radial basis function partition of unity methods for pricing vanilla basket options[END_REF]. Also, an other method for pricing American options is the operator splitting method introduced by Ikonen and Toivanen in [START_REF] Ikonen | Operator splitting methods for option pricing[END_REF] for American put option under the Black-Scholes model, and the method is studied by Toivanen in [START_REF] Toivanen | Numerical valuation of european and american options under kou's jumpdiffusion model[END_REF] under the Kou model. In this study, we consider penalty method for solving (LCP) problem ( 8), so we use penalty term defined by

δC u(τ, x) + δ -q(x)
where 0 < δ ≪ 1 is a small regularization parameter, C ≥ rK is a positive constant and q(x) = K -Ke x then we replace problem (8) by the nonlinear PIDE

∂ τ u(τ, x) -Pu(τ, x) - δC u(τ,x)+δ-q(x) = 0, u(0, x) = g(x), (10)
for (τ, x) ∈ R × (0, T ] and also we need to impose boundary conditions [START_REF] Buhmann | Radial basis functions[END_REF].

Approximation based on RBF-PUM

To avoid mesh generation, in recent years meshless techniques have attracted the attention of researchers. In a meshless method a set of scattered nodes is used instead of meshing the domain of the problem. The technique of RBFs is one of the most recently developed meshless methods that has attracted attention of many researchers in recent years. Radial basis functions are an efficient tool for solving multivariate scattered data interpolation problems. To interpolate an unknown

function f ∈ C(Ω) whose values on a set X = {x 1 , x 2 , ..., x N } ⊂ Ω ⊂ R d are known, a function of the form s f,X (x) = N j=1 λ j φ( x -x j ), (11)
is chosen, where φ : Ω × Ω → R is a radial basis function and • is the Euclidean norm. In Table 1 some globally supported RBFs are listed which are commonly employed in the literature. The positive constant c appearing in RBFs is called the shape parameter which dictates the flatness of the radial basis function and also has a key role on the convergence rate of the approximations and the condition number of the coefficient matrices. For more details about basic properties and types of radial basis functions, compactly and globally supported and also their wide applications in scattered date interpolations, the interested reader would be referred to recent works in this topic [25,40,[START_REF] Buhmann | Radial basis functions[END_REF][START_REF] Buhmann | Radial basis functions: theory and implementations[END_REF].

Name of function Definition Thin plate splines (TPS) (-1) k+1 r 2k log(r) Gaussian (GA) exp(-c 2 r 2 ) Invers multiquadrics (IMQ)

1 √ (cr) 2 +1 multiquadrics (MQ) (cr) 2 + 1 Conical splines r 2k+1
Table 1: Some well-known functions that generate globally supported RBFs RBFs methods have many advantages but generally coefficients matrix are dense, hence complexity of computations is expensive. In this section we introduces partition of unity method (PUM) which puts a given set of local approximation spaces together to produce a conforming global approximation.

Let Ω ⊂ R d be a bounded set, and let a covering {Ω j } M j=1 of the region Ω such that Ω ⊂ M j=1 Ω j . Also, we define

∀x ∈ Ω I(x) := {j | x ∈ Ω j }, card(I(x)) ≤ K,
where the constant K is independent of the number of patches M .

Definition 1. Let Ω ⊂ R d be a bounded set. Let {Ω j } M
j=1 be an open and bounded covering of Ω. This means all Ω j are open and bounded and Ω is contained in their union. Set δ j = diam(Ω j ) = sup x,y∈Ω j x-y 2 . A family of non-negative functions {w j (x)} M j=1 with w j (x) ∈ C k (R d ) is a k-stable partition of unity respect to cover of Ω j if:

1. supp(w j ) ⊆ Ω j 2. j w j (x) = 1 for x ∈ Ω 3. For every α ∈ N d 0 with |α| ≤ k there exist a constant C α > 0 such that for all 1 ≤ j ≤ M D α w j L∞(Ω j ) ≤ C α δ α j .
The weight functions w j are constructed by using Shepard's method [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF] as follow:

w j (x) = φ j (x) k∈I(x) φ k (x) , j = 1, 2, ..., M (12)
where φ j (x) are compactly supported functions with support on Ω j . To guarantee non-negativity and compact support in Ω j , we define in (12)

φ j (x) = φ x -c j r j , j = 1, 2, . . . , M, (13)
where {c j } M j=1 and {r j } M j=1 are the centers and radiuses of the circular, spherical or hyper-spherical patches {Ω j } M j=1 and where φ is one of the compactly supported functions with minimal degree described in [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions minimal degree[END_REF]Corololary 9.14]. Here we consider the Wendland compact supported function [START_REF] Wendland | Scattered Data Approximation[END_REF] φ(r) = (1r) 4 + (4r + 1), which belongs to C 2 for the construction of the weight functions.

The global approximation function s u (x), with x ∈ Ω, to the function u(x) is constructed as

s u (x) = M j=1 w j (x)s u,j (x) = j∈I(x) w j (x)s u,j (x),
where s u,j are local interpolants such that s u,j (x i ) = u(x i ) for each node x i ∈ Ω j . Then, the global PU approximant inherits the interpolation property, i.e. s u (x i ) = u(x i ). Using the cardinal basis functions

ψ j (x i ) = 1, i = j 0, i = j i = 1, 2, . . . , N, (14)
the local interpolant s u,j (x) is an RBF-approximant of type s u (x) = Ψ(x)u, on Ω j where Ψ(x) = [ψ 1 (x), ψ 2 (x), • • • , ψ N (x)].
For time dependent function u(t, x), we construct the global approximant built up from local RBF interpolants of type

s u (t, x) = N j=1 ψ j (x)u j (t) = Φ(x)A -1 u(t). (15)
For j ∈ {1, . . . , M }, let J(Ω j ) := {k | x k ∈ Ω j } be the set of indices of the node points that belong to the patch Ω j . For such patch Ω j , the local RBF approximation is given by

s u,j (t, x) = k∈J(Ω j ) ψ k (x)u k (t)
where s u,j (t, x k ) = u k (t) for all nodes x k ∈ Ω j and ψ k are cardinal basis functions. Hence, in the RBF-PUM, we obtain the global approximant for the time-dependent function u(t, x)

s u (t, x) = j∈I(x) w j (x)s u,j (t, x) = j∈I(x) k∈J(Ω j ) w j (x)ψ k (x)u k (t). ( 16
)
Now, for deriving an error estimation for partition of unity method, we need the following definitions Definition 2. For a bounded region Ω fill distance on X = {x 1 , x 2 , ..., x N } define as following

h X,Ω = sup x∈Ω min x j ∈X x -x j 2 ,
which can be explained as largest distance that for every x ∈ Ω there is one point x j with this distance. Definition 3. A subdomain Ω j ⊆ R d satisfies an interior cone condition if there exists an angle θ ∈ (0, π 2 ) and a radius γ > 0 such that, for all x ∈ Ω j , a unit vector ξ(x) exists such that the cone

C(x, ξ(x), θ, γ) = {x + λy : y ∈ R d , y 2 = 1, y T ξ(x) ≥ cos(θ), λ ∈ [0, γ]} is contained in Ω j . Definition 4. Suppose Ω ⊆ R d is bounded and X = {x 1 , x 2 , ..., x N } ⊆ R d are
given. An open and bounded covering {Ω j } M j=1 is called regular for (X, Ω) if the following properties are satisfied 1. For each x ∈ Ω, the number of subdomains Ω j , with x ∈ Ω j is bounded by a global constant C, 2. Each subdomain Ω j satisfies an interior cone condition, 3. The local fill distances h X N j ,Ω j are uniformly bounded by the global fill distance h X,Ω .

For each positive definite function φ ∈ C k (R d ) and each area Ω ⊂ R d there is a function space N φ (Ω), the native Hilbert space [START_REF] Schaback | A unified theory of radial basis functions[END_REF]. The smoothness of φ inherited from the native space via

N φ (Ω) ⊆ C [ k 2 ]
(Ω). Also, for getting the full approximation order, the weak form of Holder continuity idea is used, and we define space C k v (R d ) as space of all functions such that their derivatives of order k satisfy

D α u(x) = O( x v
2 ) for x 2 → 0. By using the above definitions the following convergence theorem is derived from [START_REF] Wendland | Scattered Data Approximation[END_REF].

Theorem 2. Suppose φ ∈ C k v (R d
) is conditionally positive definite of order m, also, let Ω ⊆ R d be open and bounded and X = {x 1 , x 2 , ..., x N } ⊆ Ω. Let {Ω j } M j=1 be a regular covering for (Ω, X) and let {w j } M j=1 be k-stable for {Ω j } M j=1 . Then the error between u ∈ N φ (Ω) and its partition of unity interpolant is bounded by

|D α u(x) -D α s u (x)| ≤ C 1 h k+v 2 -|α| X,Ω |u| N φ (Ω)
for all x ∈ Ω and all |α| ≤ k 2 .

Discretization

In this section, we split the PIDE operator P into two parts as follows

Pu(x, τ ) = Du(x, τ ) + J u(x, τ ),
where D is a differential operator and J is an integral operator defined by

Du(τ, x) = 1 2 σ2 ∂ 2 x u(τ, x) + (r -q - σ2 2 -λ(ǫ)ξ(ǫ))∂ x u(τ, x) -(r + λ(ǫ))u(τ, x) (17)
J u(τ, x) = |y|≥ǫ u(τ, x + y)ν(y)(dy).

For the spatial variable x, the domain is restricted to a bounded set Ω = [x min , x max ] and generate uniform grid points on the truncated region. For given positive integers N and M , let

△τ = T M and △x = xmax-x min N -1
. Now, let τ m = m△τ for m = 0, 1, ..., M and x n = x min + (n -1)△x for n = 1, ..., N . We define u(x n , τ m ) to be the price at the nth spatial and mth temporal node.

To use the RBF-PUM approximation for spatial discretization of PIDE problems ( 4) and ( 10), we need to compute the effect of applying a spatial differential operator D at the interior node points. Let L be a derivative operator of order α. By applying of operator L in RBF-PUM approximation for interior point we get

Lu(τ, x) = ∂ |α| ∂x α u(τ, x) = i∈I(x) k∈J(Ω i ) ∂ |α| ∂x α [w i (x)ψ k (x)] u k (τ ) = i∈I(x) k∈J(Ω i )   β≤α α β ∂ |α-β| ∂x α-β w i (x) ∂ |β| ∂x β ψ k (x)   u k (τ ). ( 19
)
For composite linear operators like D, we sum up the contributions from each term, and denote the global differentiation matrix under operator D by D.

To numerical approximation of the integral operator J at each grid point (τ k , x i ), we first use transformation z = x + y, then we have

J u(τ k , x i ) = |z-x i |≥ǫ u(τ k , z)ν(z -x i )dz.
Now, let N as a constant to cut off the domain R to a bounded interval (x -N , x N ) where the values x -N and x N will be chosen based on standard financial arguments, and zero belongs to [x -N , x N ], and split the integral operator to the following four integrals

J u(τ k , x i ) = x l i x -N u(τ k , z)ν(z -x i )dz + x i -ǫ x l i u(τ k , z)ν(z -x i )dz + x L i +1 x i +ǫ u(τ k , z)ν(z -x i )dz + x N x L i +1 u(τ k , z)f (z -x i )dz (20)
where l i and L i are the greatest integers such that

x l i < x i -ǫ < x l i +1 , x L i < x i + ǫ < x L i +1 .
Now these integrals can be approximated using the trapezoidal rule and for evaluating the second and third integrals, we need the values of u(τ k , x iǫ) and u(τ k , x i + ǫ) to apply the trapezoidal rule. For doing this since

x l i < x i -ǫ < x l i +1 and x L i < x i + ǫ < x L i +1 , the values u(τ k , x i -ǫ)
and u(τ k , x i + ǫ) can be approximated by using the linear interpolation on the intervals [x l i , x l i +1 ] and [x L i , x L i +1 ] respectively. Finally, by using approximate values for integrals, the integral matrix under integral operator J is derived and denoted by J.

To evaluate an European option, let

U m := [u(τ m , x 1 ), u(τ m , x 2 ), • • • , u(τ m , x N )]
⊤ (here ⊤ means transpose) as an approximate value of the solution which can be obtained by using the following Crank-Nicolson, Leap-Frog (CNLF) time stepping scheme

U m+1 -U m-1 2△τ = D( U m+1 + U m-1 2 ) + JU m , 1 ≤ m ≤ M -1 ( 21
)
where D is the differentiation matrix associated with the differential operator and J is integral matrix corresponding to the integral operator. For m = 0 we apply the following scheme

U m+1 -U m △τ = D( U m+1 + U m 2 ) + JU m . (22)
Also, the initial vector is defined by

U 0 = [g(x 1 ), g(x 2 ), • • • , g(x N )] ⊤ where g(•)
is the payoff function (6). One of the sources of error can arise when we use initial function (6) at grid points, since this function is discontinuous in their first derivatives. A useful notion in the implementation of numerical method is that the value of a function on a grid represents average value of the function over the surrounding grids rather than its value sampled at each grid point by

g(x i ) ≈ 1 △x x i + △x 2 x i -△x 2 g(x)dx,
and this makes the payoff function smooth at the strike price K, and we use this technique to improve the accuracy of RBF-PU numerical method especially near the strike price. For American option, after spatial discretization of differential and integral operators of ( 10), the solution can be obtained by using the following (CNLF) time stepping scheme

U m+1 -U m-1 2△τ = D( U m+1 + U m-1 2 ) + JU m + F m , 0 ≤ m ≤ M -1 (23)
where the initial vector is

U 0 = [g(x 1 ), g(x 2 ), • • • , g(x N )] ⊤ where g(•)
is the payoff function ( 6), and

F m = [ δC u(τ m , x 1 ) + δ -q(x 1 ) , δC u(τ m , x 2 ) + δ -q(x 2 ) , • • • , δC u(τ m , x N ) + δ -q(x N ) ] ⊤ .

Numerical Results

In this section, we present numerical results to show efficiency and performance of RBF-PU method and discuss more about convergence of proposed technique. All computations are carried in MAT-LAB with a 3.6 GHz Corei3 processor. Although the proposed scheme works for all radial basis functions from Table 1 we choose the multiquadric radial basis function for all the numerical experiments. We present absolute errors to check the efficiency and accuracy of RBF-PU method and for evaluating absolute error we have used reference solutions used in literature. Also, to show the practical performance of proposed numerical method, we report condition number of linear systems and CPU time in second.

The accuracy of RBF-methods where the differentiation matrix contains the differentiations of the shape function highly depends upon the shape parameter c of the basis functions, which is responsible for the flatness of the functions. For smooth problems, the best accuracy is typically achieved when c is small, but then the condition number of the linear system becomes very large. Figure 1 displays the dependence of root mean square (RMS) error on the size of the shape parameter for European and American options for some common test examples. For computing RMS error, we selected N = 2048 and M = 2048 as numbers of space and time steps and evaluated option price as a reference value for both European and American options in the region of interest

[ K 2 , 3K 2
]. We only considered parameter sets 2, 6 and 7 and other sets make similar results. Figure 1 shows that RBF-PU method is less sensitive to the change of shape parameter and we can derive that c = 2 is a good choice for European options and c = 3 is a good choice for American options.

European options

This section of numerical results consists in computing European put options with different sets of parameter values used in literature [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF][START_REF] Wang | Robust numerical valuation of European and American options under the CGMY process[END_REF][START_REF] Almendral | Accurate evaluation of European and American options under the CGMY process[END_REF] given in Table 2 in order to verify the convergence and stability of RBF-PU method combined with three time level (CNLF) discretization. We compare the numerical values with the Carr-Madan formula [START_REF] Carr | Option valuation using the fast fourier transform[END_REF] as the reference option price.

We employed parameter set 1 to evaluate European put option and numerical results have been reported in Table 3. In [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF] an implicit method is used to evaluate European and American options under Lévy model. The proposed methods solve partial integro-differential equations for European options and linear complementarity problems for American options via an operator splitting method involving linear systems with tridiagonal matrices. Also, in [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF] parameter set 1 has been used for European put option and results reported in Table 3. Comparing results in Table 3 with Table 3 of [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF] confirm that RBF-PU method is more accurate than proposed method in [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF] and needs less time and space steps to reach an specific error level. Also, CPU times and condition numbers presented in Table 3 show the efficiency of RBF-PU method and overcome ill conditioning issue of global RBF methods. Parameter set 1 has been used in [START_REF] Wang | Robust numerical valuation of European and American options under the CGMY process[END_REF] and the resulting (PIDE) is then solved using a preconditioned BiCGSTAB method coupled with a fast Fourier transform. Results reported in Table 3 compare favorably with those in Table 8 of [START_REF] Wang | Robust numerical valuation of European and American options under the CGMY process[END_REF]. [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF][START_REF] Wang | Robust numerical valuation of European and American options under the CGMY process[END_REF] 0.1 0 0 1 8.8 9.2 1.8 0.25 10 CGMY, parameter set 2, see [START_REF] Almendral | Accurate evaluation of European and American options under the CGMY process[END_REF] 0.1 0 0 1 5 5 0.1 1 1 G M Y T K CGMY, parameter set 3, see [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF] 0.07 0.2 0 0.5 7 9 1.4 0.25 10 VG, parameter set 4, see [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF] 0.06 0 0 1 8 6 0 1 10 CGMY, parameter set 5, see [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF] 0.06 0 0 1 10 8 0.6 1 10 CGMY, parameter set 6, see [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF][START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF] 0.04 0.15 0 0.5 10 10 1.2 0.25 10 VG, parameter set 7, see [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF] 0.06 0 0 4 8 8 0 0.5 10

r σ q C G M Y T K CGMY, parameter set 1, see
Table 5: Parameter setup for American put options.

We considered parameter set 6 given in Table 5 for American put option and results presented in Table 7. Also, in [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF] parameter set 6 has been used for American put option and results reported in Table 5. Comparing results confirm that RBF-PU method needs less time and space steps to reach an specific error level, and definitely it will be faster than the proposed numerical scheme discussed in [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF].

To show the convergence of the time discretization, we numerically investigate the behavior of global temporal errors as a function of ∆τ which is defined by

Maximum Error = max j |u m (T, x j ) -u 3 000 (T, x j )|, (24)
where u m is the numerical solution of European and American options at the spatial nodes x j ∈ (log( 1 2 ), log( 3 2 )) associated with S j ∈ ( K 2 , 3K 2 ) after m time steps, and u 3 000 the corresponding solution for m = 3 000 used as an approximation for the exact solution. Figure 2 displays the global temporal errors versus ∆τ sequence of seven increasing time steps m, namely 8, 16, 32, 64, 128, 256, 512 and N = 512 spatial nodes for American and European options. As a first main observation, the global temporal errors decrease monotonically as M increases or equivalently ∆τ decreases. Concerning the actual convergence behavior, it is easy to see from Figure 2 that the temporal errors as a function of ∆τ are bounded from above in each set of parameter by C(∆τ ) p with some moderate constants C where p ≈ 1 for more cases.

Also, to validate that the RBF-PU method is numerically well-behaved, we present the maximum error convergence with respect to the number of grid points over the region ( K 2 , 3K 2 in Figure 3. The maximum errors presented in 3 show that errors decrease monotonically as N increases or equivalently ∆x decreases. Concerning the actual convergence behavior, it is easy to see from Fig- 3 that the temporal errors as a function of N bounded from above in each set of parameter by C(N ) -p with some moderate constants C where p ≈ 2 for more cases.

Binary barrier and butterfly options

One of the advantages of the RBF-PU method as presented in this study is that one can deal with complicated payoffs such as binary barrier and butterfly options without affecting the accuracy. The numerical experiments confirm this fact. Therefore, in this subsection, we consider a European binary down-and-out barrier option that is used to pricing a credit default swap by Cariboni and Schoutens in [START_REF] Cariboni | Pricing credit default swaps under Lévy models[END_REF]. The payoff function of the binary down-and-out barrier option is given by [START_REF] Lord | A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes[END_REF] Table 6: Price of American put options with parameter sets 3, 4 and 5 with truncated domain is [-3, 3] and ǫ = 0.09. method discussed in this paper is rather applicable and provides accurate numerical results for non-convex payoffs. A butterfly spread is the result of buying call options with strikes K 1 and K 2 respectively, and selling two call options with strike price

g(x) = 1 {x>ln( B K )} ,
K 3 = K 1 +K 2 2
, Then, the payoff function will be g(x) = max(e x -K 1 , 0) -2 max(e x -K 3 , 0) + max(e x -K 2 , 0).

A European style butterfly is appropriate when the investor thinks large moves the price of the underlying asset are unlikely. It usually has a middle strike K 3 around the spot price K of the underlying asset and it pays off when the price of the underlying asset stays close to K. Numerical results for European and American butterfly options with parameters:

r = 0.1, σ = 0, q = 0, K = 1, T = 1, C = 1, G = 3, M = 5, Y = 0.5
obtained by RBF-PU method are given in Table 9 including option price and absolute error at S = K = 1, CPU times in second and condition number of arisen linear system of equations. Also, in Table 10, we report l ∞ -errors to show the accuracy of RBF-PU method in the region of interest [ K 2 , 3K 2 ] not only on the specific point such as S = K, and for evaluating these errors we use Carr-Madan formula [START_REF] Carr | Option valuation using the fast fourier transform[END_REF]. The option value for European butterfly is plotted in Figure 4.

Greeks

In this section our aim is to evaluate the efficiency of the RBF-PU method for computing the Greeks ∆ = ∂v ∂S and Γ = ∂ 2 v ∂S 2 which are important for hedging the option. Delta measures the sensitivity of the option price to a fluctuation in the underlying asset price while Gamma is the speed at which the Delta changes with respect to this movement. We plotted in Figures 6 and7 the Delta and of the number of patches in RBF-PU method also leads to an improvement of the approximation. The effect of the time discretization is measured by studying the temporal error. For the American and European option cases where the CNLF scheme is combined with penalty method, we conclude that the rate of convergence is of at least first order in time for more cases. Also, Figure 3 confirms that our scheme is stable and second order convergent in space. The experiments also demonstrated that European and American option prices with error around 1.0e -4 can be computed in than one second on a PC. Thus, the developed method is very fast and accurate.

Figure 1 :

 1 Figure 1: Root mean square (RMS) error as function of the shape parameter c for European and American options.

Figure 2 :

 2 Figure 2: Temporal error vs. ∆τ with N = 512 spatial nodes for American and European options.

Figure 3 :

 3 Figure 3: Maximum error convergence with respect to the number of points for American and European options. The Maximum error is measured over the region ( K 2 , 3K 2 ) with N = 3 000 and M = 512 as reference solution.

Figure 4 :

 4 Figure 4: Left: The option value for European butterfly option. Right: The option value for European binary down-and-out barrier option.

Figure 5 :

 5 Figure 5: Left: The option value for American butterfly option. Right: Early exercise boundaries for American butterfly option with Y 0.1 and Y = 0.7.

Figure 6 :

 6 Figure 6: Delta and Gamma of European butterfly option against the reference solutions.

Figure 7 :

 7 Figure 7: Delta and Gamma of American butterfly option against reference solutions.

Figure 8 :

 8 Figure8: The uniform grid with 1024 computational nodes. The partition of unity is performed using 256 circular patches of equal size. The RBF-PU coefficient matrix that was obtained on this grid. It contains only 8166 non-zero elements.

Table 2 :

 2 Parameter setup for European put options.

Also we used parameter set 2 for different values of Y to evaluate European put option and r σ q C

  where the value B is a barrier level and is the indicator function.

			S = 10, parameter set 3
	N	M	Price	Error	time(s) cond
	128	64	1.216663	5.2e-2	0.006	10.11
	256	128	1.251458	1.8e-2	0.04	24.46
	512	256	1.268218	1.2e-3	0.26	50.71
	1024	512	1.269537	2.1e-5	1.64	99.56
	Ref. Price [32] 1.269516
			S = 10, parameter set 4
	N	M	Price	Error	time(s) cond
	128	64	0.489189	6.6e-3	0.006	2.33
	256	128	0.497642	1.7e-3	0.03	2.71
	512	256	0.496501	6.2e-4	0.26	3.70
	1024	512	0.495835	3.9e-5	1.67	6.91
	Ref. Price [32] 0.495875
			S = 10, parameter set 5
		M	Price	Error	time(s) cond
	128	64	0.812406	6.9e-2	0.006	4.11
	256	128	0.875805	6.1e-3	0.04	9.07
	512	256	0.881678	2.7e-4	0.25
	1024	512	0.881959	1.0e-5	1.63	38.57
	Ref. Price

Table 8 :

 8 Price of European binary down-and-out barrier option with truncated domain is [-2, 2] and ǫ = 0.07.

			S = 100
	N	M	Price	Error time(s) cond
	128	64	0.945296 1.1e-2	0.006	7.12
	256	128	0.943196 9.7e-3	0.04	15.30
	512	256	0.943335 9.8e-3	0.25	25.51
	1024	512	0.933494 2.4e-5	1.63	50.43
	Ref. Price [31] 0.933474

Table 11 :

 11 Condition number of localized RBF-PU method and global RBF method for different sets of parameter and N = 1024 as space step size.

		Condition number of RBF-PU Condition number of global RBF
	parameter set 1	617.86	2.97e20
	parameter set 2, Y = 0.1	12.57	1.31e20
	parameter set 2, Y = 0.5	81.38	1.27e20
	parameter set 2, Y = 0.8	259.40	7.26e20
	parameter set 3	99.56	1.44e21
	parameter set 4	6.91	1.15e20
	parameter set 5	38.57	6.46e20
	parameter set 6	45.70
	parameter set 7	12.88	1.49e21

Acknowledgements

The authors would like to thank professor Ariel Almendral Vazquez, Senior Research Scientist SAND, Norwegian Computing Center, for the provided reference solution for European and American butterfly options and some fruitful comments.

 [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF]4.671721 4.392112 4.139657 Price in [START_REF] Wang | Robust numerical valuation of European and American options under the CGMY process[END_REF] 4.3716708

Table 3: Price of European put options with parameter set 1 and truncated domain is [-5, 5] and ǫ = 0.07.

numerical results have been reported in Table 4. Parameter set 2 has been used in [START_REF] Almendral | Accurate evaluation of European and American options under the CGMY process[END_REF] and in this paper partial integro differential equation is discretized in space by the collocation method and in time by an explicit backward differentiation formula for evaluation of European option. Numerical results reported in Table 2 of [START_REF] Almendral | Accurate evaluation of European and American options under the CGMY process[END_REF], and by comparison accuracy and efficiency of RBF-PU method are confirmed.

American options

In this section, we consider parameter sets given in Table 5 for evaluation of American put option. For all test cases, the reference prices for the American options are calculated using the CONV method described in [START_REF] Lord | A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes[END_REF] and given in [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF]. For American put option, we consider three different sets 3, 4 and 5 which in these cases we let Y ∈ (1, 2), Y = 0 and Y ∈ (0, 1) respectively, and numerical results are given in Table 6 including option price and absolute error at S = 10, CPU times in second and condition number of linear system of equations. Also, to increase the rate of convergence in time and space for American option problems, we implement a Richardson extrapolation method defined by

where u( △τ 2 , △x 2 ) refers to the approximate solution by using the RBF-PU method with time step size △τ 2 and spatial step size △x 2 and u(△τ, △x) refers to the approximate solution with time step size and spatial step size △x.

In [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF] exponential time integration is used for solving American option pricing problem under infinite activity Lévy models with parameter sets 3, 4 and 5 and the performance of this scheme is compared with the Crank-Nicolson method and an implicit-explicit method in conjunction with an extrapolation. In [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF] the exponential time integration scheme is superior over the entire range of errors and reaches an accuracy of the order of 1.2e -4 in 22.197 seconds, whereas the Crank-Nicolson algorithm takes 58.669 seconds and the implicit-explicit extrapolation scheme reaches the same accuracy in 94.297 seconds. By comparing CPU times reported in Table 6, it is clear that RBF-PU method due to the sparsity of linear system is more faster than the proposed methods in [START_REF] Nisha | Numerical pricing of American options under infinite activity Lévy processes[END_REF]. The corresponding parameters under the variance gamma model are given by r = 0.0421, σ = 0, q = 0, C = 1 υ , B = 50, T = 1, K = 100

, Y = 0 where σ = 0.20722, υ = 0.50215, and θ = -0.22898. These parameters are also taken in [START_REF] Cariboni | Pricing credit default swaps under Lévy models[END_REF]. We use the truncated domain [-2, 2] of the log price.

In [START_REF] Cariboni | Pricing credit default swaps under Lévy models[END_REF] authors used the numerical method developed in [START_REF] Hirsa | Pricing American options under variance gamma[END_REF] for pricing a binary down-andin barrier option under the variance gamma model. The price of the binary down-and-in barrier in [START_REF] Cariboni | Pricing credit default swaps under Lévy models[END_REF] is 0.0253 at S = 100. Then the reference value of the binary down-and-out barrier option is 0.933474 at S = 100 according to the relationship between two options [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF].

results for European binary down-and-out barrier option obtained by RBF-PU method are given in Table 8 including option price and absolute error at S = K = 100, CPU times in second and condition number of arisen linear system of equations. These parameters have been used in [START_REF] Lee | Tridiagonal implicit method to evaluate European and Amerioptions under infinite activity Lévy models[END_REF] for evaluating values of European binary down-and-out barrier option by the implicit method with three time levels and results presented in Table 6. Comparing results clearly confirm that for specific error level, the proposed RBF-PU method needs less number of time and spatial steps. Also, the option value for European binary down-and-out barrier option is plotted Figure 4. Now, we consider European and American butterfly options to check the efficiency and accuracy of developed meshfree method. The reason for including a butterfly option is to show that the ] and ǫ = 0.15 for European option and ǫ = 0.09 for American option.

Gamma functions against the reference solutions for European and American butterfly options. In these figures reference solutions are derived from Carr-Madan formula [START_REF] Carr | Option valuation using the fast fourier transform[END_REF].

Properties of the RBF-PU discretization matrices

If we use global RBF approximation coefficient matrix is dense ordinary and condition is very Using of unity method of RBF introduces sparsity in the discretization matrices. In Figure 8 the sparsity pattern of the discretized matrix in RBF-PU is shown. Also, we considered global RBF method for discretization of (4) and condition number of discretized matrix is given in Table 11 for different sets of parameter and space step sizes. Comparing results with condition number of localized RBF-PU method shows that linear system of equations arisen in RBF-PU is well conditioned.

Conclusion

We proposed RBF-PU method for spatial discretization of PIDE to price American and European options under Lévy model. Also, penalty method is used to solve the free boundary problem arisen in pricing American option. Then, CNLF is applied for time discretisations and we treat the local term using an implicit step and the nonlocal term using an explicit step to avoid to avoid the inversion of the nonsparse matrix. These result in a linear algebraic system with a sparse matrix that has a small condition number. The shape parameter in the RBF affects the accuracy and stability of the numerical methods. Figure 1 confirms that RBF-PU method is less sensitive to the change of shape parameter. An increase of the number of nodal points and correspondingly an increase


Luca Pesaresi
email: luca.pesaresi@imperial.ac.uk

Jason Armand

Christoph W Schwingshackl

Loïc Salles

Chian W Wong

O Tatin

An advanced underplatform damper modelling approach based on a microslip contact model

Keywords: Underplatform damper -nonlinear dynamics -microslip model -friction nonlinearity Fx, Fy, Fz Contact force coefficients X, Ỹ, Z Contact displacement coefficients F nl, P Nonlinear forces, external forces. Q Harmonic displacements M.C.K Mass, damping and stiffness matrices p Pressure distribution q Tangential relative displacement

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

One of the most common causes of failure for gas turbines is high-cycle fatigue (HCF) driven by large resonant stresses [ ]. Underplatform dampers are devices used by aeroengine manufacturers to reduce the risk of HCF failure, since they dissipate the vibrational energy while at the same time acting as a seal between blades to increase the thermodynamic efficiency [ , ]. Due to the nonlinear nature of the contact forces, characterised by stick-slip and separation, the dynamic behaviour of the blades becomes very nonlinear making its analysis much more complex. Continuous research has been conducted in this field for the development of new numerical damper models which a empted to improve the fidelity and/or computational cost [ , , , ]. Most of the modelling approaches proposed rely on frequency domain solvers based on the multi-harmonic balance method in conjunction with continuation techniques, as they represents an efficient approach when only the steady-state response is of interest [ , , ]. Various contact models have been used to simulate the interaction between the platform and the blades, which can be classified as macroslip models [ , , , ] or microslip models [ , , ]. Initially, the macroslip contact models were only D [ ], and subsequently they were extended to allow the relative normal motion between the damper and platform [ , , ]. In [ ], a fully coupled D macroslip contact model, capable of simulating a coupled in-plane D motion and the normal relative motion was developed. With regards to the microslip models, various concepts have been developed [ , , ] with the most common being a parallel elasto-plastic springs arrangement, initially proposed by Iwan in [ ] for D tangential motion, and then recently extended to allow relative normal motion in [ , ]. When a microslip contact model is employed, normally only a single element describes the behaviour of the whole contact interface and specifically-designed experimental set-ups are used to tune the element [ , , ]. is approach works well for non-conforming contact surfaces (cylindrical dampers), since the contact is normally very localised. However, it could lead to an over-simplification for conforming contacts (wedge dampers) since the kinematics of the contact, which is partially lost using one single microslip element, plays a significant role on the nonlinear dynamic behaviour [ ]. A different approach to model microslip, which allows an improved description of the contact kinematics, is to discretise the contact area with a grid of macroslip contact elements [ , , , ]. When this approach is chosen, a significant microslip behaviour can be observed only when a pressure gradient is present at the contact interface as shown in [ ]. However, in real contacts, frictional dissipation can also occur with a nominally uniform pressure distribution, since on the scale of the asperities no contact is truly conforming due to the surface roughness [ ]. For this reason, the overall microslip behaviour observed could be seen as the combination of two microslip processes with different scales, a larger scale microslip due to non-uniform pressure distribution, and the microscale microslip at the asperity level. In this study a new modelling approach is presented, which is characterised by the discretisation of the contact area with a grid of contact elements, which are tuned to reproduce the microscale energy dissipation. For this scope, a new D microslip contact model based on a modified Valanis elasto-plastic model [ , ] was used, and it was then tuned numerically using a semi-analytical contact solver [ ]. e modelling approach was then evaluated against an underplatform damper test rig recently developed by the authors [ , ], showing some promising capabilities for the prediction of the nonlinear dynamics of the blades.

. MODELLING APPROACH

e modelling approach proposed in this paper is based on the previous work presented in [ , ], but introduces a new contact element and a empts to increase the fidelity of the contact description at the platform-damper interface. is proposed approach is well suited to model wedge dampers, as well as other structures which have friction nonlinearities caused by conforming contacts. e damper is included in the model explicitly, using an FE description which allows to account for its flexibility and inertia properties. is aspect is important for dampers which are relatively flexible and that operate close to a "stuck" condition, since neglecting the local flexibility might lead to inaccuracies in the calculation of the contact forces.

e damper and the blades FE models are then coupled using a series of contact elements newly developed and de-scribed in details in the next section. To allow an accurate description of the contact properties, the damper-platform interface is discretised with a grid of contact elements, each of those representing the contact behaviour of a portion of the contact as shown in Fig. . Due to the manufacturing tolerances, even nominally conforming surfaces can be characterised by a lack of conformity which, together with non-uniform normal loading might lead to a non-uniform pressure distribution as shown in the example of Fig. . is variation of pressure is taken into account by explicitly defining a different initial pre-load on each element depending on its position on the contact interface, and on the area of the associated patch. e microslip behaviour on the microscale level, caused by the surface roughness, as well as the contact stiffness associated with each portion of the contact, is taken into account by tuning the parameters of each contact element, as described in section .

e numerical analysis is based on an exisiting multiharmonic balance (MHBM) solver with continuation coupled with a model reduction technique, included in the code FORSE (FOrced Response SuitE) and discussed in detail in [ , , ]. Due to the friction forces arising at the contact interface, the equations of motion for the blades and damper are nonlinear, and can be wri en in the following form:

Mx(t) + Cx(t) + Kx(t) + F nl (x(t), x(t)) = P(t). ( )
where M, C, K are the mass, damping, and stiffness matrices. P are the external excitation forces and F nl are the nonlinear contact forces dependent on the relative motion of the interacting nodes at the interface. According to the MHBM, the response for each DOF of the system can be expressed as a Fourier series truncated at the n th harmonic:

q(t) = Q 0 + n j=1 Q c j cos m j ωt + Q s j sin m j ωt ( )
where Q 0 , Q c j and Q s j are the harmonic coefficients for each DOFs. When equation is inserted in equation , and the harmonic terms balanced with a Galerkin projection, a system of equations to determine all harmonic components can be obtained in the frequency domain as follows:

Z(ω)Q + F(Q) -P = 0 ( ) where Q = Q 0 , Q c 1 , Q s 1 , . . . , Q c n , Q s n
T is a vector of harmonic coefficients, F(Q) is the vector of nonlinear forces and P is the external excitation. Z(ω) is the dynamic stiffness matrix of the system. To reduce the size of the problem, a model reduction based on a FRF matrix representation is applied [ ], and the harmonic coefficients are then calculated with a Newton-Raphson iterative solver.

.

D MICROSLIP CONTACT ELEMENT

A new contact element, which enables to simulate D contact motion while taking into account the microslip effects due e Valanis model was then succesfully employed in [ ] to model the frictional hysteretic behaviour of a lap joint subject to harmonic excitation. In the present study, the model was modified to allow the inclusion of a variable time-dependent normal load F z (t) as originally proposed in [ ]. With this approach, the friction limit is time-dependent allowing a full coupling between the tangential and normal relative displacements at the contact.

erefore, the tangential friction force is defined as follows:

F = E 0 q 1 + λ E 0 q |q | (E t q -F) 1 + k λ E 0 q |q | (E t q -F) = f (q, q, F, λ) ( )
where E 0 corresponds to the tangential contact stiffness of the initial loading, E t is the slope of the macroslip region, k is a parameter which characterises the microslip transition, q is the contact tangential displacement, and λ is a parameter defined as:

λ = E 0 µF z (t) 1 -k E t E 0 ( )
One of the advantage of the Valanis formulation is that, by defining a positive value for E t , it permits to reproduce contact hysteresis loops where the macroslip region has a slope. is contact behaviour has been observed experimentally in friction joints that are subject to fre ing wear, as reported in [ , ].

To simulate a full D motion, two Valanis elements were used along the x and y direction of the local coordinate system defined for the element "i" as shown in do describe a D in-plane motion even if there is no coupling between the two directions. For each contact element "i", a set of parameters needs to be defined: E 0 , E t , k, the friction coefficient µ, the normal contact stiffness k n , and the initial element pre-load N 0 = σ i A i , with σ i and A i being the local pressure and contact area associated with the element. In this study, a simplification was made assuming the parameters E 0 and k n constant with time, but the model is flexible to allow the inclusion of different load-dependent stiffness laws in the future.

An alternate frequency-time method (AFT) [ ] is used to allow an easier calculation of the contact forces which are nonlinear with displacements, using a separate routine in the time-domain. A simplified flow-chart for the contact forces calculation is shown in Fig. .  e contact relative displacements X, Ỹ and Z of the mating nodes between the platform and the damper are expressed as harmonic coefficients in the local coordinate system. An inverse Fourier transformation is then applied, and the non-X, Ỹ, Z

x(t), y(t), z(t) linear contact forces are then evaluated, using an implicit Euler's method for the tangential forces. e contact forces for the three directions in the local coordinate system are then transformed back to the frequency-domain, Fx , Fy , Fz and this force calculation is repeated until convergence is reached for the Newton-Raphson scheme. e parameter "k", which controls the level of "microslip" in the contact, can vary from "0" (infinite microslip) to "1" (no microslip), and the effect of this variation under a constant normal load is shown in Fig. .  One of the drawbacks of using the Valanis model for friction contacts, is that unlike in other joint models, the initial loading slope, which corresponds to E 0 (point A in Fig. ), is different from the "re-loading slope" (point B in Fig. ). In fact, from equation , for E t = 0, the re-loading slope is equal to 2E 0 /(1 + k), and therefore it is always higher than E 0 apart from the limit case k = 1, which is not defined. Since only the steady-state solution is of interest here, the initial load is always neglected, and for a correct analysis, the dependence of the contact stiffness on k must be taken into account when specifying the input parameters.

F z = N 0 + k n z, if N 0 + k n z > 0 F z = 0, if N 0 + k n z < 0 F x = 0, if F z = 0 F y = 0, if F z = 0 F x = F valanis (x, x, F x , F z ), if F z > 0 F y = F valanis (y, y, F y , F z ), if F z > 0 F x (t), F y (t), F z (t) Fx , Fy , Fz
To evaluate the capability of the new contact model under coupling conditions, a tangential harmonic motion coupled to an in-phase normal motion with a separation event is shown in Fig. . As a result of this coupling, the tangential friction force is not constant anymore at the sliding limit, but follows the curve µN(t). When the amplitude of the .

Contact displacement [mm] ×10 -3 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Contact forces [N]

TEST CASE -UPD RIG

To validate the proposed modelling approach, an underplatform damper (UPD) rig recently developed at Imperial College London (see [ , ]) was chosen as the main test case. e UPD rig is an experimental set-up that allows the investigation of the effect of UPDs on blade-like structures under a controlled lab environment. e assembly of the rig can be seen in Fig. .  Two pseudo beam-like blades are fixed to a common base, simulating a rigid disk, and are clamped via a hydraulic cylinder to a large inertia block.

e damper is a wedge type, [ ], which has a triangular cross section with a characteristic angle. Unlike a real high pressure turbine blade, the aerofoil is replaced by a straight rectangular cross-section beam, but still maintaining vibration modes similar to a real blade. e centrifugal load is simulated via a pulley system with calibrated masses, and the excitation is provided by an electrodynamic shaker (Data Physics V ) a ached near the root of the blade to minimise its impact. A non-contact measurement system, which employs two single point laser Doppler vibrometers (Polytec OFV-) was used to capture the dynamic response of both blades near the tip (see Fig. ). A stepped sine test was performed in a narrow frequency band around the mode of interest, always keeping the excitation force constant within a tolerance range to control the nonlinearity of the system. Different tests were performed with increasing excitation forces in order to progressively activate the friction nonlinearity at the damper-platform interface and characterise the dynamics of the system.

.

CONTACT ELEMENT TUNING

One of the fundamental steps of the proposed modelling approach is the tuning of the contact parameters for the newly-developed friction contact element. In fact, including more accurate information about the contact nonlinearity, which is the driver for the nonlinear behaviour of the system, is very important to improve the prediction capabilities of the proposed modelling approach. e approach proposed here for the contact element tuning is mainly based on a semi-analytical contact solver described in section . , but fully experimental approaches based on specifically-designed test set-ups could be used as well. e advantage of using a semi-analytical contact solver relies on the possibility to simulate any rough surface, without any assumption on the statistical distribution of the asperities. In addition, the relatively high speed of computation allows to overcome the problems of FE rough contact simulations, which would require prohibitively large models to account for a good representation of the surface features.

. Semi-analytical contact solver based on the Boundary Element Method (BEM)

A BEM-based contact solver, previously developed and validated by the authors [ , ], is used here to perform a refined analysis of the damper-platform friction interface and tune the new contact element. e contact solver uses the projected conjugate gradient method [ ] and a discreteconvolution fast Fourier transform to accelerate the computation. It assumes the elastic half-space body description, which makes it possible to use the Boussinesq and Cerruti potentials [ , ] to compute the surface elastic deflections in the normal and tangential directions from the pressures and shear tractions in the contact area. Equation gives the component of normal displacement u z due to a pressure distribution p. Equation is the discretised form of Eq. on a regular grid of N x × N y points. Similar equations are used to compute the the tangential displacements. e first step of the algorithm is to solve the normal contact problem using the conjugate gradient method, a er which the tangential problem can be solved using the Coulomb friction law to bound the shear distribution in the slipping region.

u z (x, y) = 1 -ν 2 πE +∞ ∫ -∞ +∞ ∫ -∞ p(ξ, η) (ξ -x) 2 + (η -y) 2 dξdη ( )
where E and ν are the Young's modulus and Poisson ratio of the material, respectively.

u z (i, j) = K zz ⊗ p = N x k=1 N y l=1 p(k, l)K zz (i -k, j -l) ( )
where ⊗ denotes the discrete convolution product and K zz (i, j) are the discrete influence coefficients [ ] that give the normal displacement resulting from unit pressure on the element centred on the grid point (i, j).

.

Pressure distribution tuning

Initially, a D non-contact topography-measuring interferometer Polytec TMS-was used to scan both contact surfaces of the damper used, revealing slightly curved profiles (see [ ]). However, due to the limited dimensions of the scanning area of the interferometer, it was not possible to measure the surfaces of the platforms. erefore, this limited the possibility to use the semi-analytical contact solver to obtain the contact pressure distribution, and an experimental approach was preferred instead. Fujifilm prescale pressure films were used at the damper-platform interface, and the prints obtained were then scanned and digitised in MATLAB. A 4 th order polynomial fit was used to approximate the pressure film prints (see Fig. ), which was then interpolated to obtain the input pre-loads for each contact element. . Microscale tuning e semi-analytical contact solver described in section . was used to numerically tune the newly-developed contact element, allowing the model to reproduce the microslip effects due to the surface roughness, as well as the characteristic stiffness of the contact. Initially, a stylus profilometer was used for both the damper and the platforms to evaluate the roughness characteristics of the contact surfaces. A patch of 250µm by 250µm was generated in MATLAB, as shown in Fig. ,using a Gaussian distribution for the asperities which was able to reproduce the experimental roughness (with R q = 0.177µm and a correlation length Cl x = 6.6µm). A contact mesh of by elements was used for the contact simulations, as it allows a good discretisation of the surface asperities, while keeping an acceptable computational time. Figure .

Variation of the parameter "k n " with contact pressure.

For the tangential problem, contact simulations were run imposing a sinusoidal displacement, and the corresponding contact hysteresis loops were obtained for different pressure levels ( . MPa, . MPa, MPa, MPa, MPa, MPa) in the range of interest, as shown in Fig. . Since both E 0 and k influence the stiffness of the contact and its microslip transition, their selection is not trivial, and an optimisation was run in MATLAB based on the interior-point algorithm described in [ ]. A set of parameters E 0 and k could be found for each pressure level, which allowed the model to perfectly capture the numerical hysteresis loops. By interpolating the plots of E 0 and k versus pressure (see Fig. and Fig. ), it was possible to specify a set of optimised parameters for each contact element at the damper-platform depending on the element pressure. A monotonic increasing trend was obtained for E 0 , which indicates the tendency of the initial loading stiffness to increase at higher pressure, whereas a more complex behaviour characterised by a minimum close to MPa was obtained for the microslip parameter k. .

RESULTS AND DISCUSSION

In this section, the prediction capabilities of the newly-proposed modelling approach are evaluated against the experimental data of the UPD rig ). In addition, similarly to the previous model, it was decided to use the pressure distribution interpolated from the pressure films, as described in section . . erefore, the difference obtained with the previous model will allow to evaluate the contribution of the newly-developed contact element, as well as the tuning procedure described in section . for a higher fidelity description of the contact properties. In fact, in the previous model [ ], the contact stiffness was assumed the same (60000N/mm 3 ) for all the contact elements, based on previous measurements performed on the D friction rig described in [ ].

e mode investigated for the present study is the first flexural out-of-phase (OOP) mode shown in Fig. , as the damper operates close to a microslip condition, whereas the in-phase mode is dominated by a large separation at the contact interface (see [ ]), and it is therefore much less interesting for the evaluation of the new approach. Some initial experiments were performed for the first flexural mode without the damper, in order to evaluate the background damping caused by the material and set-up configuration including the shaker. A loss factor of η = 0.04% was measured and it was included as modal damping in the linear model. e nonlinear measurements with the damper in place were performed using a constant damper load of N, which allowed a good conformity at the contact interface, and an increasing shaker excitation force from . N up to N (shaker limit). e nonlinear FRFs of the experiments, new modelling approach, and old modelling approach are all compared in Fig. , and the excitation and response location can be seen in Fig. .  anks to an improved description of the stiffness distribution at the contact interface, the resonance frequency of the nearly-linear excitation case for the new model is slightly closer to the experiments, being . % higher compared to the . % of the previous model. Both models tend to overestimate the resonance frequency of the experimental OOP mode, which lies at . Hz. A possible explanation for this behaviour could be an overestimation of the contact area in the model caused by the uncertainty in the polynomial approximation of the experimental pressure films. With regards to the amplitude estimation, the new model significantly improves the prediction accuracy compared to the previous model. In having an amplitude percentage difference of -. %, + . %, + . %, + . %, + . %, + . % for the increasing excitation levels compared to the + . %, + . %, + . %, + %, + . % and + . % of the previous model. erefore, taking into account the microscale miscroslip dissipation due to the surface roughness with a new microslip element lead to a significantly improved amplitude estimation, highlighting the promising capabilities of the new model for accurate dynamic predictions. e slight overestimation of the experimental nonlinear amplitudes observed is still acceptable, since factors such as the additional damping which might have been introduced by the damper pulling wire were not included in model. More investigations need to be performed for an improved evaluation of the damper-platform contact area, since this is one of the main factors for an accurate prediction of the modes resonance frequencies when the damper is in place.

. CONCLUSIONS

In this study, a new modelling approach for underplatform dampers was presented and evaluated against the experimental data of a recently-developed test rig. e new model employs a dense discretisation of the damper-platform interface with a grid of D microslip contact elements. For this scope, a new D contact element based on a modified Valanis model, was developed and preliminary tests highlighted its capabilities for representing strong contact nonlinearities. e key idea of the model proposed, is to capture both the contribution of a macroscale microslip, due to non-uniform pressure distribution, as well as the microscale microslip due to the surface roughness. A tuning procedure for the microslip elements is also proposed, based on the numerical simulations of the normal and tangential behaviour of a representative contact patch having similar roughness characteristics to the real contact. e new modelling approach proved to be very effective in capturing the amplitude levels of the nonlinear FRFs, leading to a significant improvement compared to the previous models based on the use of macroslip contact elements. ese findings suggest the relevance of the microscale dissipation, as well as of the improved description of the contact stiffness, in the overall frictional dissipation of the damper. In addition, the modelling approach proposed, only relies on the measurement of some contact properties, such as the friction coefficient and roughness, and therefore it represents a promising tool for the prediction of the nonlinear dynamic behaviour of blades constrained by friction dampers.

for giving permission to publish it.

is work is part of a collaborative R&T project SILOET II P . which is co-funded by Innovate UK and Rolls-Royce plc and carried out by Rolls-Royce plc and the Vibration UTC at Imperial College London.

Figure .

 . Figure . Discretisation of the damper-platform contact interface.

FigFigure .

 . Figure . New D microslip element schematic.

Figure .

 . Figure .Flow-chart of the contact element routine.

Figure .

 . Figure .Contact hysteresis loops with varying microslip parameter "k".

Figure .

 . Figure . Contact hysteresis loops with increasing normal displacement.

Figure .

 . Figure . Underplatform damper test rig lab set-up.

Figure .

 . Figure .Pressure film digitised and approximated with a polynomial fit.

Figure .

 . Figure . Simulated rough contact patch.

Figure .

 . Figure . Simulated contact hysteresis loops at different pressure levels. e friction Coulomb limit was set to be µ = 0.6, since this was the average value measured for the same material in previous experimental tests with a D friction rig developed at Imperial College London [ ]. e hysteresis loops of Fig. were then used to tune the contact element input parameters E 0 and k.Since both E 0 and k influence the stiffness of the contact and its microslip transition, their selection is not trivial, and an optimisation was run in MATLAB based on the interior-point algorithm described in [ ]. A set of parameters E 0 and k could be found for each pressure level, which allowed the model to perfectly capture the numerical hysteresis loops. By interpolating the plots of E 0 and k versus pressure (seeFig. and Fig. ), it was possible to specify a set of optimised parameters for each contact element at the damper-platform depending on the element pressure. A monotonic increasing trend was obtained for E 0 , which indicates the tendency of the initial loading stiffness to increase at higher pressure, whereas a more complex behaviour characterised by a minimum close to MPa was obtained for the microslip parameter k.

Figure . Figure .

 .. Figure .Variation of the parameter "E 0 " with contact pressure.

  , and are comapared to the previous model used in [ ]. To allow a comparison with the previous results, the finite element model of the blades and damper were kept the same (see Fig. ) as in [ ], as well as the subset of contact nodes used to discretise the damper-platform interface (see Fig.

Figure .

 . Figure . First flexural OOP mode of the blades and damper.

Figure .

 . Figure . Contact friction elements at the damper-platform interface.

Figure .

 . Figure .Nonlinear FRFs at increasing excitation levels.

Figure .

 . Figure .Comparison of peak nonlinear amplitudes at increasing excitation levels.

  L. Pesaresi, L. Salles, A. Jones, J.S. Green, and C.W. Schwingshackl. Modelling the nonlinear behaviour of an underplatform damper test rig for turbine applications. Mechanical Systems and Signal Processing, : -, . [ ] C. M. Firrone and S. Zucca. Modelling friction contacts in structural dynamics and its application to turbine bladed disks. In J.Awrejcewicz, Numerical Analysis -eory and Application, INTECH , pages -, . [ ] C. Putignano, M. Ciavarella, and J.R. Barber. Frictional energy dissipation in contact of nominally flat rough surfaces under harmonically varying loads. Journal of the Mechanics and Physics of Solids, ( ): -, . [ ] K. C. Valanis. Fundamental consequences of a new intrinsic time measure. plasticity as a limit of the endochronic theory. Arch. Mech., : -, . [ ] L. Gaul and J. Lenz. Nonlinear dynamics of structures assembled by bolted joints. Acta Mechanica, ( ): -, . [ ] J. Armand, L. Salles, and C. W. Schwingshackl. Numerical simulation of partial slip contact using a semianalytical method. In ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, pages V T A -V T A . American Society of Mechanical Engineers, . [ ] L. Pesaresi, L. Salles, R. Elliot, A. Jones, J. S. Green, and C. W. Schwingshackl. Numerical and experimental investigation of an underplatform damper test rig. Applied Mechanics and Materials, : -, . [ ] E. P. Petrov. A High-Accuracy Model Reduction for Analysis of Nonlinear Vibrations in Structures With Contact Interfaces. Journal of Engineering for Gas Turbines and C. W. Schwingshackl, E. P. Petrov, and D. J. Ewins. Measured and estimated friction interface parameters in a nonlinear dynamic analysis. Mechanical Systems and Signal Processing, : -, April . [ ] J. Armand, L. Pesaresi, L. Salles, and C.W. Schwingshackl. A Multiscale Approach for Nonlinear Dynamic Response Predictions With Fre ing Wear. Journal of Engineering for Gas Turbines and Power, ( ), . [ ] I. A. Polonsky and L. M. Keer. Fast methods for solving rough contact problems: A comparative study. J. Tribol., ( ): , . [ ] K. L. Johnson. Contact Mechanics. Cambridge University Press, . Cambridge Books Online. [ ] A. E. H. Love. A Treatise on the Mathematical eory of Elasticity. A Treatise on the Mathematical eory of Elasticity. at the University Press, . [ ] C .W. Schwingshackl. Measurement of Friction Contact Parameters for Nonlinear Dynamic Analysis. In Proceedings of the Society for Experimental Mechanics IMAC Conference, . ] R. H. Byrd, J. C. Gilbert, and J. Nocedal. A trust region method based on interior point techniques for nonlinear programming. Mathematical Programming, ( ): -, Nov .

					Dampers in Turbomachinery Blading : Part -Stick-
					Slip Contact Kinematics. Journal of Engineering for Gas
					Turbines and Power,	(April	): -,	.
					[ ] B.D. Yang and C.H. Menq. Characterization of d con-
					tact kinematics and prediction of resonant response of
					structures having d frictional constraint. Journal of
	[ ] B.A.Cowles. High cycle fatigue in aircra gas turbine -an industry prospective. International Journal of Fracture,	Sound and Vibration, [ ] E. P. Petrov and D. J. Ewins. Analytical Formulation of ( ): -, .
	: -,	.		Friction Interface Elements for Analysis of Nonlinear
	[ ] J. H. Griffin. A review of friction damping of turbine blade vibration. International Journal of Turbo and Jet Engines, : -, .	Multi-Harmonic Vibrations of Bladed Disks. Journal of Turbomachinery, ( ): , . [ ] C. H. Menq, J. Bielak, and J.H. Griffin. e influence of
	[ ] K. Y. Sanliturk, D. J. Ewins, and A. B. Stanbridge. Un-derplatform Dampers for Turbine Blades: eoretical Modeling, Analysis, and Comparison With Experimen-tal Data. Journal of Engineering for Gas Turbines and	microslip on vibratory response, part I: A new microslip model. Journal of Sound and Vibration, ( ): -, June . [ ] G. Csaba. Forced Response Analysis in Time and Fre-
	Power,	( ): ,	.	quency Domains of a Tuned Bladed Disk With Friction
	[ ] E. P. Petrov. Explicit Finite Element Models of Fric-tion Dampers in Forced Response Analysis of Bladed Disks. Journal of Engineering for Gas Turbines and Power,	Dampers. Journal of Sound and Vibration, , July . [ ] K.Y. Sanliturk and D. J. Ewins. Modelling Two-( ): -
	( ):	,	.		Dimensional Friction Contact And Its Application Using
	[ ] C. M. Firrone, S. Zucca, and M. M. Gola. underplatform dampers on the forced response of bladed e effect of disks by a coupled static/dynamic harmonic balance	Harmonic Balance Method. Journal of Sound and Vibra-tion, : -, . [ ] W. D. Iwan. A distributed-element model for hysteresis
	method. International Journal of Non-Linear Mechanics,	and its steady-state dynamic response. Journal of Applied
	( ): -,	.		Mechanics, ( ): -,	.
	[ ] E. Cigeroglu, N. An, and C. H. Menq. Forced Response Prediction of Constrained and Unconstrained Structures Coupled rough Frictional Contacts. Journal of Engi-neering for Gas Turbines and Power, ( ): , . [ ] L. Panning, W. Sextro, and K. Popp. Spatial Dynamics of Tuned and Mistuned Bladed Disks with Cylindrical and Wedge-Shaped Friction Dampers. International Journal of Rotating Machinery, : -, . [ ] M. Krack, L. Salles, and F. ouverez. Vibration predic-tion of bladed disks coupled by friction joints. Archives of Computational Methods in Engineering, pages -, . [ ] M. Krack, L. Panning, and J. Wallaschek. A high-order harmonic balance method for systems with distinct states. Journal of Sound and Vibration, ( ): -, . [ ] L. Salles, B. Staples, N. Hoffmann, and C. W. Schwing-shackl. Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions. Nonlinear Dynamics, ( ): -, ( ): , . [ ] F. Pourahmadian, H. Ahmadian, and H. Jalali. Modeling and identification of frictional forces at a contact inter-face experiencing micro-vibro-impacts. Journal of Sound and Vibration, ( ): -, . [ ] Power, [ ] M. Rajaei and H. Ahmadian. Development of gener-alized Iwan model to simulate frictional contacts with variable normal loads. Applied Mathematical Modelling, ( ): -, . [ ] C. Gastaldi and M. M. Gola. An Improved Microslip Model for Variable Normal Loads, pages -. Springer International Publishing, Cham, . [ ] G. Csaba. Modelling of a microslip friction damper sub-jected to translation and rotation. In Proceedings of ASME Gas Turbine and Aeroengine Congress and Exhibition, vol-ume , page V T A , . [ ] C. Gastaldi and M. M. Gola. On the relevance of a mi-croslip contact model for under-platform dampers. In-ternational Journal of Mechanical Sciences, -: -, . [ ] L. Pesaresi, M. Stender, V. Ruffini, and C. W. Schwing-shackl. DIC Measurement of the Kinematics of a Friction Damper for Turbine Applications, pages -. Springer International Publishing, Cham, . [ ]
	.
	[ ] J. H. Griffin. Friction Damping of Resonant Stresses in
	Gas Turbine Engine Airfoils. Journal of Engineering for
	Power,	: -,	.
	[ ] B. D. Yang and C. H. Menq. Characterization of Contact
	Kinematics and Application to the Design of Wedge

[

ACKNOWLEDGMENTS

e authors are grateful to Innovate UK and Rolls-Royce plc for providing the financial support for this work and