1,721,137 research outputs found
Credit Portfolio Losses
Full text linkApplied ScienceElectrical Engineering, Mathematics and Computer Scienc
Pricing multi-asset options with sparse grids
No full textMulti-asset options are based on more than one underlying asset, in contrast to standard vanilla options. A very significant problem within the pricing techniques for multi-asset options is the curse of dimensionality. This curse of dimensionality is the exponential growth of the complexity of the problem when the dimensionality increases, because the number of unknowns to solve simultaneously grows exponentially. Modern computer systems cannot handle this huge amount of data. In order to handle the multi-dimensional option pricing problem, the curse of dimensionality has to be dealt with. The sparse grid solution technique is one of the key techniques to do this. The sparse grid technique divides the original problem into many smaller sized sub-problems, which can be handled efficiently on a modern computer system. Because every sub-problem is independent of all others, this technique is parallelisable at a high efficiency rate. This means, that every sub-problem can be solved simultaneously. However, because of the dimensionality, the size of the sub-problems may remain too large to solve and should be parallelised further. The main restriction to the application of the sparse grid method is that the mixed derivative of the solution of a multi-dimensional option pricing problem has to be bounded. Because of the typical non-differentiability of the final condition of the option pricing problem, this restrictions has to be taken seriously. In the first part of this thesis, it is shown, experimentally, that indeed the sparse grid technique does not lead to a satisfactory accuracy without the use of advanced techniques. If a coordinate transformation is used, the accuracy increases significantly. This transformation aligns the non-differentiability along a grid line. Coordinate transformations are not applicable to any type of multi-asset option, which seriously restricts the sparse grid solution technique for real life financial applications. Sometimes, however, it is not necessary to use it, because the non-differentiability is already aligned with grid line. These types of options are the options based on the best or worst performing underlying asset. The boundary conditions of these contracts are unknown en henceforth these options are computed and analysed with a second alternative method in this thesis. This method arises from the risk-neutral expectation valuation of the final condition which can be written as a multi-dimensional integral over the transition density. By use of a discrete Fourier transform, we can solve this integral efficiently. The fast Fourier transform is a fast algorithm to compute the discrete Fourier transform. This algorithm serves as the basis for a sophisticated algorithm to parallelise the computation of the discrete Fourier transform, by dividing the transform in several parts. In this thesis, a complete parallel algorithm which does not require communication between the sub-problems is developed, which subdivides the problem in a sophisticated way. In combination with the sparse grid technique, the numerical results have a satisfactory accuracy.Electrical Engineering, Mathematics and Computer Scienc
Efficient Pricing of Early: Exercise and Exotic Options Based on Fourier Cosine Expansions
No full textIn the financial world, two tasks are of prime importance: model calibration and portfolio hedging. For both tasks, efficient option pricing is necessary, particularly for the calibration where many options with different strike prices and different maturities need to be priced at the same time. Therefore, a fast yet accurate pricing method is a necessity for banks and trading companies. Nowadays three groups of pricing methods are being used in the financial industry and academia, that is, Monte--Carlo methods, partial (integro-)differential equation (PIDE) methods, and numerical integration methods, where the option price is modeled as the discounted expected value of the payoff at maturity. The latter type of methods is attractive from both practice and research point of view, as the fast computational speed, especially for plain vanilla options, makes it useful for calibration at financial institutions. Usually numerical integration techniques are combined with the Fast Fourier transform or Hilbert transform, and therefore, the numerical integration methods are often referred to as the `transform methods'. Representatives of transform methods are the Carr--Madan method (Carr, Madan, 1999), the CONV method (Lord et.al. 2008) and the Hilbert transform method (Feng, Linetsky, 2008). A recent contribution to the transform method category is the COS method proposed in Fang, Oosterlee (2008, 2009), that is, an option pricing method based on the Fourier cosine expansions. It departs from a truncated risk--neutral formula, in which the conditional density function is recovered in terms of its characteristic function, by Fourier cosine expansions. This method can be used for asset processes as long as the characteristic function of the conditional density function is known, or can be approximated. For processes where the density function and its derivatives are continuous functions with respect to the underlying asset, the COS method exhibits an exponential convergence rate. Our research work is based on the COS method, which has been used for vanilla European option pricing (Fang Oosterlee, 2008), vanilla early--exercise option pricing and barrier option pricing (Fang, Oosterlee, 2009). The motivation of this thesis is to further improve the robustness of the COS method, make it efficient for non--Levy models, and extend it to different types of exotic options. The point of departure of this thesis is to improve the robustness of the COS method for call option pricing with early-exercise features, as presented in Chapter 1, where the call option prices are obtained from put option prices, in combination with the put--call parity and put--call duality relations, which are incorporated into our pricing algorithm at each early--exercise date to recover the Fourier coefficients and to compute the continuation value. The robustness of the pricing methods is demonstrated by error analysis, as well as by a series of numerical examples. In Chapter 2, the acceleration of option pricing by the COS method on the Graphics Processing Unit (GPU) is presented. After a brief discussion of the GPU and its potential for option pricing, we will study different ways of GPU implementation, followed by three examples of GPU acceleration, the so-called multiple strike option pricing, option pricing under hybrid models where the characteristic function is derived from a Riccati ODE system, and the example of Bermudan option pricing. Influence of data transfer between host and device is also discussed in this chapter. Extension of COS method to early--exercise option pricing with an Ornstein--Uhlenbeck (OU) process is explained in Chapter 3. OU processes for commodity derivatives, either with or without seasonality functions, are non--Levy processes and more computationally expensive within the COS framework, as compared to Levy processes. First of all, an accurate pricing algorithm is given, which can be used for all OU processes with different types of seasonality functions. Then, based on a detailed error analysis, a more efficient pricing method is proposed, which reduces the computing time from seconds to milliseconds. However, this new method is not advocated for all parameter settings. The conditions under which the basis point accuracy can be ensured is derived by error analysis. In the numerical part, the accuracy and efficiency of these two pricing methods are compared, and the conditions we derived from error analysis are further verified by several numerical experiments. In Chapter 4, we present an efficient pricing method for American--style swing options, based on Fourier cosine expansions. Here we assume that the holder of the swing option has the right, but not the obligation, to buy or sell a certain amount of commodity, such as gas and electricity, at any time before the expiry of the option, and more than once. Moreover, a recovery time is added between two consecutive exercises in which exercise is not allowed. Our pricing method is based on the Bellman principle, leading to a backward recursion procedure in which the optimal exercise regions are determined at each time step, after which the Fourier coefficients can be recovered recursively. Our method performs well for different underlying processes, different swing contracts and different types of recovery time. The pricing methods for European and early--exercise Asian options (ASCOS) are shown respectively in Chapters 5 and 6. In Chapter 5, we present an efficient option pricing method for Asian options written on different types of averaged asset prices, but without early--exercise features. In our method, the characteristic function of the average asset is recursively recovered, with the help of Fourier expansions and Clenshaw--Curtis quadrature. Then it is used in the risk--neutral formula to get the Asian option price. Exponential convergence rate is observed for most Levy processes, which is also supported by a detailed error analysis. Advantages of our pricing algorithm are that as the number of monitoring dates increases, the method stays robust and the computing time does not increase significantly, as shown in the numerical results. Our pricing method for early--exercise Asian options is presented in Chapter 6. In this case, the Fourier cosine coefficients of the option price are recursively recovered by Fourier transform and Clenshaw--Curtis quadrature. Then these coefficients are inserted into the risk--neutral formula, which, in the early--exercise Asian case, is a two--dimensional integration, to get the option value. The chain rule from probability theory is also needed in our algorithm to factorize the joint conditional density functions. An exponential convergence rate in the option price, as derived in a detailed error analysis, is observed from various numerical experiments. Factors of approximately hundred of speedup are achieved on the GPU. Conclusions and insight into future research are to be found in Chapter 7. In this thesis, efficient pricing methods for different early--exercise and exotic options, based on the Fourier cosine expansions, are presented, followed by an error analysis and numerical results, from which we see that the COS method is an efficient, robust and flexible method for pricing different types of option products, for different asset models, and is suitable for GPU acceleration. It is a promising tool for financial calibration and dynamic hedging in practice.Applied mathematicsElectrical Engineering, Mathematics and Computer Scienc
Equity and Foreign Exchange Hybrid Models for Pricing Long-Maturity Financial Derivatives
No full textModelling derivative products in Finance usually starts with the specification of a system of Stochastic Differential Equations (SDEs), that corresponds to state variables like stock, interest rate, Foreign Exchange (FX) rate and volatility. By correlating the SDEs for the different asset classes one can define the hybrid models, and use them for pricing multi-asset derivatives. Even if each of the individual SDEs yields a closed-form solution, a non-zero correlation structure between the processes may cause difficulties for efficient product pricing. Typically, a closed-form solution of hybrid models is not known, and numerical approximation by means of Monte Carlo (MC) simulation or discretization of the corresponding Partial Differential Equations (PDEs) has to be employed for pricing. The speed of pricing European derivative products is crucial, especially for the calibration of the SDEs. Several theoretically attractive SDE models, that cannot fulfil the speed requirements, are not used in practice.Numerical AnalysisElectrical Engineering, Mathematics and Computer Scienc
The COS Method: An Efficient Fourier Method for Pricing Financial Derivatives
Full text linkWhen valuing and risk-managing financial derivatives, practitioners demand fast and accurate prices and sensitivities. Aside from the pricing of non-standard exotic financial derivatives, so-called plain vanilla European options form the basis for the calibration of financial models. As any pricing and risk management system has to be able to calibrate to these plain vanilla options, it is important to be able to value these options quickly and accurately. By means of the risk-neutral valuation formula the price of any option, without early exercise features, can be written as an expectation of the discounted payoff of this option. Starting from this representation one can apply three types of numerical methods to calculate the price itself: Monte Carlo simulation, numerical solution of the corresponding partial-(integro) differential equation (P(I DE) and numerical integration. An important aspect of research in Computational Finance is to continuously improve the performance of the pricing methods. In this dissertation we present a novel and efficient option pricing method based on the integration representation for various financial derivatives. The method is called the COS method, because the key idea is to replace the probability density function, appearing in the risk neutral valuation formula, by its Fourier-cosine series expansion. Fourier-cosine series coefficients have an elegant closed-form relation with the characteristic function (i.e. the Fourier transform of the underlying density function), which is readily known, for example for exponential Lévy processes, and more general for the broader class of regular affine processes. For European options, presented in Chapter 2 of this thesis, the risk-neutral valuation formula appears as an inner product of the probability density function and the payoff function. Approximating the density function by its Fourier-cosine series expansion transforms the inner product into combinations of products of cosine basis functions and the (payoff) function which is known analytically. The method's convergence is therefore directly connected to the convergence of the Fourier-cosine series approximation of the density function. For smooth density functions, we show that the convergence is exponential. Since the computational complexity grows only linearly with the number of terms in the expansion, the COS method is optimal in error convergence and in computational complexity for European options. In other chapters of this thesis, the applicability of the COS method has been generalized to pricing options with early-exercise features and discretely-monitored barrier options, as well as to the calibration of Credit Default Swaps, under exponential Lévy processes. Furthermore, an efficient pricing method based on the COS method has been developed for pricing early-exercise options under the (two-dimensional) Heston stochastic volatility dynamics. The main insight for generalizing to the options with early-exercise features under Lévy processes, as presented in Chapter 3, is that between any two adjacent early-exercise dates the COS formula for European options can be applied, and that the series coefficients of the option values at each time lattice can be recovered recursively from those of the payoff function. This recursion can be written as a matrix-vector product with the matrix being the sum of two special matrices: a Hankel matrix and a Toeplitz matrix. Multiplication of a vector to these special matrices can be written as circular convolution of two vectors and can therefore be computed rapidly by means of the Fast Fourier Transform (FFT) algorithm. The overall computational complexity for early-exercise (and discretely monitored barrier) options with M exercise dates is O((M-1)Nlog2N), with N being the number of terms in the cosine expansion. It is then shown in Chapter 4 that the COS method can be used to calibrate Credit Default Swaps (CDSs), that are the basic building blocks of the credit risk market. The (bid-ask) spread for a CDS depends on the default probability of the underlying reference entity. In the approach adopted here, the credit default spreads are related to a series of survival/default probabilities with different maturities. These survival probabilities can be viewed as prices of binary down-and-out barrier options (without discounting). As such, the COS method for discrete barrier options can be used to recover the probabilities. A special scheme has been developed to reduce the computational time by computing several survival probabilities simultaneously. The method's capabilities have been demonstrated by the calibration to quotes of the constituents of the iTraxx Series. Finally, in Chapter 5, the COS method has been extended by quadrature rules to efficiently deal with two-dimensional pricing problems, originating for early-exercise options under the Heston stochastic volatility model. We focus especially on parameter values for which the volatility can reach zero. This is a nontrivial situation for which many other numerical schemes fail, including the finite difference PDE schemes. The problem is related to parameter values for which the Feller condition is not satisfied, and for which the pricing problem is close to singular. The variance of the stock process follows a noncentral chi-square distribution, and, for some combinations of the relevant parameters, the density function values tend to infinitely large numbers when the variance approaches zero. We handle this problem by a transformation from the variance domain to the log-variance domain. The two-dimensional discrete pricing formula is then obtained by applying the Fourier-cosine series expansion approximation in the log-stock dimension and a high-order quadrature rule in the log-variance dimension. The overall computational complexity is almost-linear in the log-stock dimension and quadratic in the log-variance dimension with very satisfactory error convergence. In this thesis we show that the COS method can price various financial derivatives under exponential Lévy and Heston stochastic volatility models. Highly efficient pricing results are presented, due to a fast error convergence and a lean computational complexity. We further determine the stability of the pricing method by a rigorous error analysis and by many numerical experiments under extreme parameter values.Applied MathematicsElectrical Engineering, Mathematics and Computer Scienc
Efficient Multigrid Methods based on Improved Coarse Grid Correction Techniques
No full textMultigrid efficiency often suffers from inadequate coarse grid correction in different prototypic situations. We select a few problems, where coarse grid correction issues arise because of anisotropic coefficients, non-equidistant or non-uniform grid stretching, or inherent indefiniteness in the partial differential equation. Most of the work in this thesis can be classified as an attempt to increase multigrid efficiency by analysing and developing novel grid coarsening techniques that ensure sufficient coarse grid correction for the multigrid algorithm. Anisotropy in discrete systems can stem from various continuous and discrete features of the problem and has to have its negative effects countered before a successful multigrid solution can be brought about. We select multidimensional stationary diffusion equation as the first important problem to be treated in this context. The work for dimensions higher than three, is aimed at developing grid coarsening strategies for discretization on rectangular hyper-grids that differ greatly in their dimensions, and thus induce the so-called grid-aligned anisotropies in the system. Coarse grids formed through standard coarsening fail to provide sufficient coarse grid correction, and alternative block relaxation techniques are expensive in high dimensions. We also investigate and test coarsening strategies with the aim that their use would allow point based relaxation to stay effective in this non-equidistant multigrid scenario. Through local Fourier analysis we also analyze w-RB Jacobi, and implement a computer program through which we compute the optimal relaxation parameters. There are three important inferences in this regard. (1) Partial (and grid dependent) coarsening strategies allow the successful use of point relaxation methods for this problem. (2) Quadrupling along a few dimensions is a very attractive partial coarsening choice. (3) Optimal relaxation parameters have a significant enhancement effect on multigrid convergence in high dimensions. The efficient solution of time-dependent multidimensional equations (discretized with implicit time-integration schemes) is also a challenge. We first use the sparse grid technique to reduce the exponential complexity of the discrete problem, and then use the d-dimensional multigrid techniques to solve the sparse grid sub problems. In this situation, i.e., with a multitude of different non equidistant grids, evaluating and using optimal relaxation parameters on the fly is not an option any more. As a multigrid solver in high dimensions depends on optimal attributes to quite a large extent, we employ the method as a preconditioner, instead of a solver. This results in a very robust and efficient multigrid preconditioned Bi-CGSTAB solver. Another coarsening strategy that we develop in this thesis is aimed at two-dimensional grids that are non-uniformly stretched. We investigated different experimental coarsening strategies. A strategy based on improving individual mesh aspect ratios of grid cells proves successful both theoretically as well as experimentally. It is based on adaptive coarsening so that on each successive coarsening step the grid cells become more square. We can also successfully use point relaxation methods with the proposed technique and get nice multigrid convergence in this case. This bit of work also has consequence for locally refined grids. Efficient multigrid techniques for the indefinite Helmholtz equation form a separate research theme included in the thesis. We employ the complex shifted operator preconditioning technique for model problems that stem from quantum mechanics applications. These model problems have strongly varying wave-numbers which perturbs the solution. The mesh size requirement in the region of this perturbation is quite demanding. This requirement can be eased by saturating the grid in that area. We find that standard coarsening in this situation works well. This is in contrast to the existing strategies where coarsening is only done in the region of refinement, until the grid is regularized. We discretize the model problems both on equidistant and also on locally refined grids, and the efficiency of the multigrid preconditioner is tested in both these situations. Experiments point to the fact that existing techniques for the indefinite Helmholtz are still not satisfactory and must be enhanced. Some conclusions and outlook mark the end of the thesis.Applied mathematicsElectrical Engineering, Mathematics and Computer Scienc
Investment decisions under uncertainties: A case of nuclear power plants
Full text linkThis thesis discusses the role of flexibility of decisions when investing in projects that are affected by economic uncertainties. It uses the theory of real options to value such investment decisions. The thesis focuses on investment decisions related to nuclear power plants, which usually are affected by several sources of economic uncertainties. It also introduces a new Monte Carlo based option pricing method, called the Stochastic Grid Bundling Method (SGBM). SGBM can be used to efficiently price exotic financial options with early exercise features, as well as real options. Investment questions, such as, under what conditions would small modular reactors would be - from an economic perspective - a better choice, when compared to large nuclear power plants; are then addressed in the thesis.Applied mathematicsElectrical Engineering, Mathematics and Computer Scienc
Fourier Methods for Multidimensional Problems and Backward SDEs in Finance and Economics
Full text linkIn this thesis we deal with processes with uncertainties, such as financial asset prices and the global temperature. We model their evolutions by so-called stochastic processes. Many of these stochastic processes are based on the Wiener process, whose increments are normally distributed. Other models may contain jump components, to model, for example, economic disasters or degradation failures. An important class of models is the Levy class, where successive increments are independent and statistically identical over different time intervals of the same length. This may give computational advantages. A well-known application of stochastic processes is in financial mathematics, where the goal is to price financial derivatives or to estimate risk measures. The underlying asset prices may be modeled by, e.g., geometric Brownian motions. More involved models, like the Variance Gamma process, are defined by jumps. In other, for instance economic, personal, or societal, contexts one may face options in the sense of real `choices'. For example, should one build a new factory now or in the future? Or should one heighten a dike today, and by how much, or in the future? These decisions are called real options and can often be related to financial options. Similar methods can be used to value them. The numerical problems we consider deal with conditional expectations. Often these problems can be connected to a partial differential equation (PDE) by a Feynman-Kac theorem. Then we can apply PDE methods, such as finite difference schemes and finite volume methods, to approximate the solutions. From the perspective of the probabilistic representation, the class of Monte Carlo methods can be beneficial for high-dimensional problems. They are based on simulated paths of the stochastic process. Besides, the expected value can often be represented by an integral, which offers opportunities to use numerical integration techniques, such as Newton-Cotes formulas. We focus on a subclass of numerical integration methods, i.e., Fourier-based methods. These `transform methods' combine a transformation to the Fourier domain with numerical integration. The probability density function of the random variables of interest is usually unknown. Its Fourier transform, i.e., the characteristic function, is however often known and can be used to approximate the corresponding density and distribution function. Our method of choice is the COS method, which is based on Fourier cosine series expansions and the characteristic function. The matrix-vector product appearing may be computed in an efficient way by using a Fast Fourier Transform (FFT) algorithm, especially when dealing with Levy processes. Besides, the use of the discrete Fourier cosine transform helps us with the approximation of the Fourier coefficients. After a general introduction of this thesis in Chapter 1, in Chapter 2 we explain the COS formula to compute conditional expectations and provide an error analysis. Then the COS method is applied to a specific class of problems: stochastic control problems, in which an agent has the possibility to affect the trend or variation of a stochastic process in such a way that his target function is maximized. For example, he can determine his consumption or savings rate. With the COS method we approximate functions by using Fourier cosine series. Similar to Fourier series they may suffer from the Gibbs phenomenon: trying to recover a function with a jump discontinuity results in undesired oscillations, even if the number of terms in the series goes to infinity. Smooth density functions give rise to a fast exponentially converging error of the COS method, but a density function with a discontinuity in one of its derivatives results in slower algebraic convergence. This is related to the Gibbs phenomenon. A remedy to improve this is by using spectral filters, which smoothen the approximations, see Chapter 3. Vanilla call and put options are based on one underlying asset. On the other hand, rainbow options are written on multiple assets and the holder may possess a `basket' of assets. The payoff of, for example, a call-on-max option, depends on the maximum of several assets. In this way, an option holder can manage his risks. Financial options on two assets are discussed in Chapter 4. Also single-asset options under the Heston model, in which an additional stochastic process for the volatility is used, are priced by the so-called 2D-COS method, developed in this chapter. Chapter 5 deals with a problem from the field of climate change economics. The future global temperature is highly uncertain and there are different damage estimates. In our model an agent can choose the consumption level of his wealth while he is subject to these uncertain climate damages. There is a trade-off between consuming now and saving for later. Economic equilibrium conditions result in a mathematical expression for the appropriate social discount rate. Here the future temperature process and the economic wealth are the uncertain processes and we combine the methods from Chapter 2 and Chapter 4 to solve the problem. Forward stochastic processes are rather well-known. Their initial value is prescribed and a prescription for the process forwards in time is given. During the last decades the backward stochastic differential equations (BSDEs) have become popular. A BSDE is stochastic differential equation for which a terminal condition, instead of an initial condition, has been specified. Its solution consists of a pair of processes, where one of the processes `steers' the other towards the terminal condition. The value of a call option can be modeled in this way as the holder of a replicating portfolio aims to end up with a certain payoff at the terminal time. Market imperfections can also be incorporated, such as different lending and borrowing rates for money, the presence of transaction costs, or short selling constraints. In Chapter 6 we extend the COS method to solve such problems and name it the BCOS method. A forward stochastic process can be approximated by different simulation schemes. The stochastic Euler scheme is a generalization of the Euler scheme for ordinary differential equations. Higher order Taylor schemes include more stochastic terms to obtain a better convergence rate. In Chapter 7 we extend our pricing and valuation methodology by using the characteristic function of these discrete schemes within the BCOS method framework. With the second-order weak Taylor scheme and a theta-time discretization we obtain second order convergence in the number of timesteps for several problems. The techniques in Chapter 7 enable us to generalize the applicability of the BCOS method to forward SDEs for which the `continuous' characteristic function is not available.Numerical AnalysisElectrical Engineering, Mathematics and Computer Scienc
Fast solvers for concentrated elastic contact problems
No full textRail transportation plays an important role in our everyday life, and there is fast development and modernization in the railway industry to meet the growing demand for swifter, safer and more comfortable trains. At the same time, the security of train operation and the maintenance of rails have to be considered. A lot of research on these issues has been carried out, among which the study of the contact between a train's wheel and the rail is particularly significant. The contact problem considers two elastic bodies. When they are pressed together, a contact area is formed where the two body surfaces coincide with each other. Moreover, an elastic field of stress, strain and displacement arises in each body. These stresses consist of normal stress (pressure) and frictional stress (traction) acting in the tangential direction. When solving the so-called {\it normal contact problem}, we search for the contact area and the pressure on it. The {\it tangential contact problem} is studied when the two bodies are brought into relative motion. If the relative velocity of the two surfaces is small, a creeping motion may be observed which is largely caused by the elastic deformation at the contact region. In those parts of the contact area where the tangential stress is small, the surfaces of the two bodies stick to each other. Otherwise, local relative sliding may occur. The research question is to find the adhesion and slip areas, and the tangential tractions. The solution methods for contact problems have been studied from the late nineteenth century, resulting in a variety of analytic and numerical approaches, w.r.t. their own specific applications. Motivated by the requirement of fast computation for involved applications such as the simulation of railway wheel-rail dynamics, we aim at developing fast numerical solvers for concentrated elastic contact problems in this thesis. Our work focuses on the contact between bodies of linear homogeneous elastic material. Moreover, it is a concentrated contact, i.e. the contact area is small compared to the dimensions of the contacting bodies. The models in use are provided by a variational formulation, which is based on a boundary element method (BEM). It gives rise to a convex optimization problem with linear or nonlinear constraints. The corresponding Karush-Kuhn-Tucker conditions provide the governing equations and contact conditions, that are numerically solved. The most time-consuming part attributes to solving a Fredholm integral of the first kind, resulting from the BEM. The corresponding Green's function expresses the relation between tractions and deformation, using a half-space approach. This integral yields linear systems with coefficient matrices that are dense, symmetric and positive definite. Moreover, they are Toeplitz in two-dimensional (2D) problems and block Toeplitz with Toeplitz blocks in three-dimensional (3D) problems. Fast computing techniques such as the fast Fourier transform (FFT) are explored. We start our work by solving the normal contact problem in Chapter 2. It is modeled by a linear complementarity problem, for which a full multigrid method (FMG) is presented. This method combines a multigrid (MG) method, an active set strategy and a nested iteration technique. It is applied to a Hertzian smooth contact and a rough surface contact. The results show the efficiency and robustness of the FMG method. Tangential contact is considered in Chapter 3 and Chapter 4. A 2D no-slip tangential problem is first studied in Chapter 3, where we mainly solve the surface integral. A fast MG method is proposed with an FFT smoother, where a Toeplitz preconditioner is constructed to approximate the inverse of the coefficient matrix. This smoother reduces many error components but enlarges some smooth error modes. Techniques such as subdomain deflation and row sum modification (RSM) are incorporated. Numerical experiments indicate rapid convergence and mesh-independence of MG with the FFT+RSM smoother. Moreover, FFT+RSM as a stand-alone solver also shows its efficiency. The complexity of these two methods is , with the number of unknowns. We work on the 3D tangential contact in Chapter 4, where a nonlinear constrained optimization problem arises. A fast solver, called TangCG, is proposed. It combines an active set strategy and a nonlinear conjugate gradient method. The most pronounced component of this method is that it employs two types of variables in the adhesion and slip areas. Techniques including the FFT and diagonal preconditioning are also incorporated. The TangCG method is tested for Cattaneo shift problems, with different amounts of slip. It dramatically reduces the computational time, compared to the state-of-art ConvexGS method. The numerical methods presented above are based on the influence coefficients (ICs) that give the relation between tractions and deformation. In Chapter 5, we investigate ICs by computing them numerically. Based on a concentrated contact setting, an elastic model is built for this purpose and a finite element method (FEM) is employed. Suggestions about the FEM meshing and element types are given, considering the accuracy and computational cost. The effects of employing the numerical ICs on contact solutions are examined. The work in this chapter provides a guidance for developing fast solvers for conformal contact problems, which typically are governed by a larger and curved contact region. With the research presented in the present PhD thesis on efficient numerical solution techniques, the numerical solution of full-scale train-rail contact problems may have come one step closer. With the research presented in present PhD thesis, and with the resulting improved numerical solution techniques, it becomes one step closer to incorporate detailed contact models in the numerical simulation of rail vehicle dynamics and in the simulation of rail and wheel wear and track deterioration.Applied mathematicsElectrical Engineering, Mathematics and Computer Scienc
Reduction of computing time for seismic applications based on the Helmholtz equation by Graphics Processing Units
No full textThe oil and gas industry makes use of computational intensive algorithms to provide an image of the subsurface. The image is obtained by sending wave energy into the subsurface and recording the signal required for a seismic wave to reflect back to the surface from the Earth interfaces that may have different physical properties. A seismic wave is usually generated by shots of known frequencies, placed close to the surface on land or close to the water surface in the sea. Returning waves are usually recorded in time by hydrophones in a marine environment or by geophones during land acquisition. The goal of seismic imaging is to transform the seismograms to a spatial image of the subsurface. Migration algorithms produce an image of the subsurface given the seismic data measured at the surface. In this thesis we focus on solving the Helmholtz equation which represents the wave propagation in the frequency domain. We can easily convert from the time domain to the frequency domain and vice-versa using the Fourier transformation. A discretization with second-order finite differences gives a sparse linear system of equations that needs to be solved for each frequency. Two- as well as three-dimensional problems are considered. Krylov subspace methods such as Bi-CGSTAB and IDR(s) have been chosen as solvers. Since the convergence of the Krylov subspace solvers deteriorates with an increasing wave number, a shifted Laplacian multigrid preconditioner is used to improve the convergence. Here, we extend the matrix-dependent multigrid method to solve complex-valued matrices in three dimensions. As the smoother, we have considered parallelizable methods such as weighted Jacobi (?-Jacobi), multi-colored Gauss-Seidel and damped multi-colored Gauss-Seidel (?-GS). The implementation of the preconditioned solver on a CPU (Central Processing Unit) is compared to an implementation on the GPU (Graphics Processing Units or graph- ics card) using CUDA (Compute Unified Device Architecture). The results show that in two dimensions the preconditioned Bi-CGSTAB method on the GPU as well as the pre- conditioned IDR(s) method on a single GPU are about 30 times faster than on a single- threaded CPU. To achieve double precision accuracy on the GPU we have used the iterative refinement in Chapter 2. However, problems of realistic size are too large to fit in the memory of one GPU. One solution for this is to use multiple GPUs. A currently widely used architecture consists of a multi-core computer connected to one or at most two GPUs. Moreover, those GPUs can have different characteristics and memory sizes. A setup with four or more identical GPUs is rather uncommon, but it would be ideal from a memory point of view. It would imply that the maximum memory is four times more than on a single GPU. How- ever GPUs are connected to a PCI bus and in some cases two GPUs share the same PCI bus, which creates data transfer limitations. To summarize, using multi-GPUs increases the total memory size but data transfer problems appear. Therefore, in Chapter 3 we consider different multi-GPU approaches and understand how data transfer affects the performance of a Krylov subspace solver with shifted Laplace multigrid preconditioner for the three-dimensional Helmholtz equation using CUDA (Compute Unified Device Architecture). Two multi-GPU approaches are considered: data parallelism and split of the algorithm. Their implementations on a multi-GPU architecture are compared to a multi-threaded CPU and single GPU implementation. The results show that the data parallel implementation suffers from communication between GPUs and the CPU, but is still a number of times faster compared to many-cores. The split of the algorithm across GPUs limits communication and delivers speedups comparable to a single GPU implementation. As a geophysical application which requires an efficient numerical method we con- sider 3-D reverse time migration with the constant-density acoustic wave equation in Chapter 4. The idea of migration in the time domain is to calculate the forward wave- field by injecting the source wavelet. Secondly, we compute the wavefield backward in time by injecting the recorded signal at the receiver locations. Subsequently, we cross- correlate the forward and backward wavefields at given timesteps. An explicit finite- difference scheme in the time domain is a common choice. However, it requires a significant amount of disk space to store the forward wavefields. The advantage of migration with a frequency domain solver is that it does not require large amounts of disk space to store the snapshots. However, a disadvantage is the memory usage of the solver. As GPUs have generally much less memory available than CPUs, this impacts the size of the problem significantly. The frequency-domain approach simplifies the correlation of the source and receiver wavefields, but requires the solution of a large sparse linear system of equations. The question is whether migration in the frequency domain can compete with a time-domain implementation when both are performed on a parallel architecture. Both methods are naturally parallel over shots, but the frequency-domain method is also parallel over frequencies. If we have a sufficiently large number of compute nodes, we can compute the result for each frequency in parallel and the required time is dominated by the number of iterations for the highest frequency. Here, GPUs are used as accelerators and not as independent compute nodes. We optimize the throughput of the latter with dynamic load balancing, asynchronous I/O and compression of snapshots. Since the frequency- domain solver employs a matrix-dependent prolongation, the coarse grid operators required more storage than available on GPUs for problems of realistic sizes. An alternative to the depth migration is least-squares migration (LSM). LSM was introduced as a bridge between full waveform inversion and migration. Like migration, LSM does not attempt to retrieve the background velocity model, however, like full wave- form inversion the modeled data should fit the observations. In Chapter 5 an efficient LSM algorithm is presented using several enhancements. Firstly, a frequency decimation approach is introduced that makes use of the redundant information present in the data. It leads to a speedup of LSM, whereas the impact on accuracy is kept minimal. Secondly, to store the sparse discretization and matrix-dependent prolongation matrices efficiently, a new matrix storage format VCRS (Very Compressed Row Storage) is presented. This format is capable of handling lossless compression. It does not only reduce the size of the stored matrix by a certain factor but also increases the efficiency of the matrix-vector computations. The study shows that the effect of lossless and lossy compression with a proper choice of the compression parameters are positive. Thirdly, we accelerate the LSM engine by GPUs. A GPU is used as an accelerator, where the data is partially transferred to a GPU to execute a set of operations, or as a replacement, where the complete data is stored in the GPU memory. We demonstrate that using GPU as a replacement leads to higher speedups and allows us to solve larger problem sizes. Summarizing the effects of each improvement, the resulting speedup can be at least an order of magnitude compared to the original LSM method.Delft Institute of Applied MathematicsElectrical Engineering, Mathematics and Computer Scienc
- …
