Bary: The Barycentric Interpolation in MatLab

Purpose

• The library "Bary" provides MatLab functions for the barycentric interpolation.

Click Here to Download the Codes


Specifications

• Name:Bary
• Author: Feng Chen
• Version: 1.0
• Language: MATLAB

Simple Example

• Equation: $$\begin{equation} \begin{aligned} &c_1(x)u''(x)+ c_2(x)u'(x)+c_3(x)u(x) = f(x), \quad \Omega = (a,b),\\ &u \text{ satisifies the general boundary condition}. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} c_1(x)=-4x^2, \quad c_2(x) = 3x, \quad c_3(x)=2. \end{equation}$$
• Exact solution and input functions: $$\begin{equation} u(x) = \sin x , \end{equation}$$ then $f(x)$ is calculated accordingly.

Quick Start

• Compiling and running:

  matlab test_bvp
  

• Output:

   
  error =
  
     1.1180e-13
  

• CPU: Intel(R) Core(TM) i5-4590S CPU @ 3.00GHz
• OS: Ubuntu 18.04.
• Release: MATLAB R2019a.

References

• Tee, T. & Trefethen, L. A Rational Spectral Collocation Method with Adaptively Transformed Chebyshev Grid Points, SIAM Journal on Scientific Computing, Volume 28, Number 5, 1798-1811, (2006).

Core Subroutines

function [P L U A] = bc_bvp_matrix(c1, c2, c3, a1, a2, b1, b2, D1, D2)

n = length(c1);
A = diag(c1)*D2 + diag(c2)*D1 + diag(c3);
A(1,:) = a1*[1 zeros(1,n-1)] + b1*D1(1,:);
A(n,:) = a2*[zeros(1,n-1) 1] + b2*D1(n,:);
[L,U,P] = lu(A); 

FFTM: Fourier Spectral Methods in MATLAB

Purpose

• The library "fftm" provides MatLab functions that are useful for solving ellipitic systems in 1D and 2D with Fourier spectral methods. It works for both complex and real data.

Click Here to Download the Codes


Specifications

• Name:fftm
• Author: Chris Chen
• Version: 2.2
• Language: MATLAB

Simple Example

• Equation: $$\begin{equation} \begin{aligned} &a(x) u - b(x) (\beta_1 u_{xx} + \beta_2 u_{yy}) = f, \quad \Omega = (0,2\pi)^2,\\ &u \text{ is periodic}. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} a(x) =-0.1+\sin x + \cos 2y, \quad b(x) = 4 + \cos x + \sin y, \quad \beta_1=2, \quad \beta_2=3. \end{equation}$$
• Exact solution and input functions: $$\begin{equation} u(x,y) = (1+\sin x + 2\cos y)(2+ \cos y), \end{equation}$$ then $f(x,y)$ is calculated accordingly.

Quick Start

• Compiling and running:

  matlab test_rfft_elliptic_solver_2dv
  

• Output:

   
  bicgstab converged at iteration 10.5 to a solution with relative residual 8.5e-11.
  
  err =
  
       5.548761450313577e-10
  

• CPU: Intel(R) Xeon(R) CPU X5550 @2.67GHz.
• OS: CentOS release 6.4 (Final).
• Release: MATLAB R2012a.

References

• Fast Fourier transform in MATLAB. Link.

Core Subroutines

%fftr_modes returns the ordering
%of the wave number when the input data in 
%the physical space are real numbers.
%High dimensional index is the tensor 
%product of 1-D indices.
%The orderings are different for odd k 
%and even k.

%k contains the ordering.
%m is the highest mode.
function [lk m] = fftr_modes(n)


lk = zeros(n,1); 

m = floor((n-1)/2);

lk(1:m+1) = [0:m]'; 

lk(n:-1:n-m+1) = -[1:m]';

Alpha-BS: A MATLAB Package of Mathematical Algorithms for the Black-Scholes Model

Purpose

• The library "alpha-bs" computes the approximation to classic and fractional Black-Scholes models for European and American options, with finite difference method, spectral method, backward Euler method, BDF2 method, Crank-Nicolson method, projected LU method, and operator splitting method.

Click Here to Download the Codes


Simple Example

• Equation: $$\begin{equation} \begin{aligned} & v_t = \frac{\sigma^2 s^2}{2} v_{ss} + (r-q) s v_s - rv, \quad (s,t) \in (0,\infty) \times (0,T],\\ & v(s,0) = g(s). \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} T =0.25, \quad r= 0.05, \quad \sigma=0.15, \quad K=50, \quad q =0, \quad S_m = 2K. \end{equation}$$
• Exact solution for European options: $$\begin{equation} \begin{aligned} & C = Se^{-q(T-t)} N(d_1) - Ke^{-r(T-t)} N(d_2), \\ & P = Ke^{-r(T-t)} N(-d_2)-Se^{-q(T-t)}N(-d_1), \end{aligned} \end{equation}$$ where $$\begin{equation} \label{x1} d_1 = \frac{\ln(S/K) +(T-t)(r-q+\frac{\sigma^2}{2})}{\sigma \sqrt{(T-t)}}, \quad d_2 = d_1 -\sigma \sqrt{(T-t)}. \end{equation}$$ In the above, $S$ is the stock price, $K$ is the strike price, $r$ is the interest rate, $q$ is the dividend rate, $T$ is the maturity span of the option, and $t$ is the spot time (e.g., $t=0$ indicates the present time). And $N(x)$ is the cumulative distribution function (cdf) of a standard normal random variable: $$\begin{equation} \label{Nx} N(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^x e^{-\frac{s^2}{2}} ds. \end{equation}$$

Quick Start

• Compiling and running:

  matlab test_bs_eu
  

• Example Output for American options:

• CPU: Intel(R) i5-4310U.
• OS: Windows Home 10.
• Release: MATLAB R2017b.

References

• A new operator splitting method for American options under fractional black–scholes models Chris Chen, Zeqi Wang, Yue Yang, Computers & Mathematics with Applications 77 (8), 2130-2144, 2019
• A new spectral element method for pricing European options under the Black–Scholes and Merton jump diffusion models, Feng Chen, Jie Shen, Haijun Yu, Journal of Scientific Computing 52 (3), 499-518, 2012
• Operator splitting methods for American option pricing, Samuli Ikonen, Jari Toivanen, Applied mathematics letters 17 (7), 809-814

Code Highlights

% Backward Euler's method for the Black-Scholes Eqaution
[v vq] = mch_bs_eu_be(pid, N, M, M, K, T, r, sigma, q, S, g, Sq);
% BDF2 method for the Black-Scholes Eqaution
[v vq] = mch_bs_eu_bdf2(pid, N, M, K, T, r, sigma, q, S, g, Sq);
% Crank-Nicolson method for the Black-Scholes Eqaution
[v vq] = mch_bs_eu_cn(pid, N, M, 4, K, T, r, sigma, q, S, g, Sq);

VERLET: Symplectic integration for Hamiltonian systems and adjoint solvers

Purpose

• The package "Verlet" solves the second-order differential equation (or a Hamiltonian system) with the velocity Verlet scheme. It is second-order accurate in both position and velocity. In addition, the velocity Verlet scheme almost preserves the Hamiltonian (oscillating around a constant mean). The package also provides the corresponding adjoint solver.

Specifications

• Name: Verlet.
• Author: Feng Chen.
• Finishing date: 07/24/2013.
• Languages: MATLAB.

Simple Example

• Equation: $$\begin{equation} \begin{aligned} & u''=\lambda u, \quad t \in (0,T],\\ & u(0) = 0, \quad u'(0) = \sqrt{-\lambda}. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} T=10, \quad \lambda=-400, \quad \delta t = 0.0001. \end{equation}$$
• Exact solution: $$\begin{equation} u(t) = \sin(\sqrt{-\lambda} t). \end{equation}$$

Quick Start

• Compiling and running:

  matlab Verlet_Driver
  

• Output:

  maximum error of position and velocity =
       5.821930990741464e-04
  

• CPU: Intel(R) Xeon(R) CPU X5550 @2.67GHz.
• OS: CentOS release 6.4 (Final).
• Release: MATLAB R2012a.

References

• Verlet integration. Link.

Code Highlights

% the velocity-Verlet time stepping.

K = -lambda;

r = u0(1);
v = u0(2);
a = -K*r;

for i = 1:nt
 rt = r; 
 r = (1-K*dt*dt/2)*r + dt*v; 
 v = (1-K*dt*dt/2)*v + dt*(K*dt*dt/4-1)*K*rt; 
end
u = [r; v];

AGATE: GPU-accelerated 3-D Chebyshev collocation methods for elliptic equations

Purpose

• The package "Agate" provides a GPU-accelerated Chebyshev collocation method for solving the single elliptic equation with Dirichlet, Neumann, or Robin boundary conditions on a 3-D rectangular domain.

Specifications

• Name: Agate.
• Author: Feng Chen.
• Finishing date: 05/30/2013.
• Languages: CUDA C++, Fortran 90.
• Required libraries: BLAS, LAPACK, CUBLAS.

Simple Example

• Equation: $$\begin{equation} \begin{aligned} & \alpha u -(\beta_1 u_{xx} + \beta_2 u_{yy} + \beta_3 u_{zz}) = f, &&\quad \text{in } \Omega, \\ & au + b\frac{\partial u}{\partial \boldsymbol{n}} = c, && \quad \text{on } \partial \Omega. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} \Omega = (-1,1)^3 , \quad \alpha =2, \quad \beta_1=3, \quad \beta_2=4, \quad \beta_3 = 5, \quad a = 1, \quad b=1. \end{equation}$$
• Exact solution and input functions: $$\begin{equation} u(x,y,z) = e^{x+y+z}, \end{equation}$$ then $f(x,y,z)$ and $c(x,y,z)$ are calculated accordingly.

Quick Start

• Compiling and running:

  cd ./Agate
  make library
  nvcc Agate_Main.cu -llibrary -llapack -lblas
  ./a.out
  

• Output:

  Nx = 2^5, Ny = 2^6, Nz = 2^7
  Device 0: Tesla M2050 of version 2.0
  intialization time: 0.0961988
  cpu time bulk: 0.110826
  cpu error bulk: 1.57341e−11
  cpu error upx: 1.53673e−11
  cpu error umx: 3.08875e−12
  cpu error upy: 1.48668e−11
  cpu error umy: 2.78527e−12
  cpu error upz: 1.2097e−11
  cpu error umz: 1.57017e−11
  cpu error du: 2.71758e−10
  gpu time HtD: 0.00597085
  gpu time bulk: 0.00152022
  gpu time BV: 0.000412896
  gpu time diff: 0.00107117
  gpu time DtH: 0.00265923
  gpu error bulk: 1.57332e−11
  gpu error upx: 1.53664e−11
  gpu error umx: 3.08864e−12
  gpu error upy: 1.48654e−11
  gpu error umy: 2.78522e−12
  gpu error upz: 1.2097e−11
  gpu error umz: 1.57021e−11
  gpu error du: 2.71349e−10
  ...
  

• CPU: Intel(R) Xeon(R) CPU X5630 @2.53GHz.
• GPU: Nvidia Tesla M2050.
• OS: CentOS release 6.4 (Final).
• Compiler: nvcc 4.2.9, gcc/gfortran 4.5.1.

References

• Feng Chen and Jie Shen. A GPU parallelized spectral method for elliptic equations in rectangular domains, Journal of Computational Physics, Volume 250, 555-564, (2013).

Code Highlight

// 
// Mars_D_Invert performs the inversion in the frequency space
// on the device.
//
// Input:
//        dd: structure for static data
//        alpha, beta: coefficient
//        d_f: right hand side
// Output:
//        d_f: solution
// Other:
//        lambda: eigenvalue
//        s: stiffness matrix
//        q: number of quadrature points
//
inline void Mars_D_Invert (Mars& dd, double alpha, 
 double betax, double betay, double betaz, double* d_f, double* wk) {
 dim3 T(8, 8, 8);  
 dim3 B((dd.qx+7)/8, (dd.qy+7)/8, (dd.qz+7)/8);
 Mars_Kernel <<<B,T>>> (alpha, betax, betay, betaz, 
  dd.qx, dd.qy, dd.qz, dd.lambdax, dd.lambday, dd.lambdaz, wk, d_f);    
}

BARIS: Fourier method for the Burgers' equation

Purpose

• The package "BariS" solves the 1-D viscous Burgers' equation with periodic boundary condition. The spatial discretization is a Fourier psuedospectral method. Three methods were proposed to calculate the nonlinear term: (i) evaluating Fourier coefficients exactly with the convolution; (ii) padding zeros in the $k$-space; and (iii) filtering with the two thirds rule. For the temporal discretization, a second-order symmetric Strang splitting was used. The linear diffusion term was integrated exactly, while the nonlinear term by a third-order Runge-Kutta scheme.

Specifications

• Name: BariS.
• Author: Feng Chen.
• Finishing date: 10/13/2012.
• Languages: MATLAB.

Simple Example

• Equation: $$\begin{equation} \begin{aligned} & u_t + uu_x = \epsilon u_{xx}, \quad x\in (0,2\pi), \\ & u \text{ is periodic}. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} T=2, \quad \epsilon=0.001, \quad N_x = 1024, \quad \delta t = 0.005. \end{equation}$$
• Initial condition: $$\begin{equation} u_0(x) = \sin(x). \end{equation}$$

Quick Start

• Compiling and running:

  matlab Bari_Driver_Serial
  

• Graphics:
  
• CPU: Intel(R) Xeon(R) CPU X5550 @2.67GHz.
• OS: CentOS release 6.4 (Final).
• Release: MATLAB R2012a.

References

• Lloyd N. Trefethen. Spectral Methods in MATLAB (Software, Environments, Tools), SIAM, (2001).
• Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, Th.A. Spectral Methods: Fundamentals in Single Domains (Scientific Computation), Springer, (2007).

Code Highlights

%%% calculate the nonlinear term with convolution 

% re-order the modes
uu = [ui(lcut+1:N); ui(1:lcut)]/N;

% convolution
uc = conv(uu,uu);

% re-order to matlab convention
m = length(uc);
ucc = [uc((m+1)/2:m); uc(1:(m-1)/2)]*2*N;

% first-order derivative
du = rfft_grad_1d(ucc, 1);

% cut off the high modes
rhs = [du(1:lcut); du(lcut+N:m)]/2;
rhs = -rhs/2;

FITZ: Parallel integration for FitzHugh-Nagumo ODEs

Purpose

• The package "Fitz" solves the FitzHugh-Nagumo equation using the parareal method. To construct serial solvers with coarse and fine resolutions, the first-order forward Euler scheme was used.

Specifications

• Name: Fitz.
• Author: Feng Chen.
• Finishing date: 07/06/2013.
• Languages: MATLAB.

Simple Example

• Equation: $$\begin{equation} \begin{aligned} & v' = v-v^3-w+I_{\text{ext}},\\ & \tau w' = v + a - bw. \end{aligned} \end{equation}$$
• Parameters: $$\begin{equation} I_{\text{ext}}=0.5, \, a=0.7,\, b=0.8, \,\tau=12.5,\, T=200, \,\Delta T =0.1, \, \Delta t = 0.05, \, \delta t=0.001. \end{equation}$$
• Initial conditions: $$\begin{equation} v(0)=1, \quad w(0)=0. \end{equation}$$

Quick Start

• Compiling and running:

  matlab Fitz_Driver_Parareal
  

• Output:

  iteration 0, 0.000000e+00, 7.040824e-01 
  iteration 1, 1.235809e-02, 6.917243e-01 
  iteration 2, 8.904816e-04, 6.908338e-01 
  iteration 3, 6.172277e-05, 6.907721e-01 
  iteration 4, 2.008283e-06, 6.907701e-01 
  iteration 5, 9.682426e-08, 6.907702e-01 
  iteration 6, 2.131734e-08, 6.907702e-01 
  iteration 7, 2.023203e-09, 6.907702e-01 
  iteration 8, 1.402405e-10, 6.907702e-01 
  iteration 9, 7.806200e-12, 6.907702e-01 
  iteration 10, 3.443912e-13, 6.907702e-01 
  sucessful return! 
  ...
  

• Graphics:
  
• CPU: Intel(R) Xeon(R) CPU X5550 @2.67GHz.
• OS: CentOS release 6.4 (Final).
• Release: MATLAB R2012a.

References

• Lions, J. L and Maday, Y. and Turinici, G. A "parareal" in time discretization of PDE's, Comptes Rendus de l'Academie des Sciences Series I Mathematics, Vol 332, (2001).
• The FitzHugh-Nagumo model. Link.

Code Highlights

% parareal iteration
for k = 2:imax+1
 % fine solver
 for i = k:n+1
  ytF = Fitz_Time(sF, yp(:,i-1), l*m); yF(:,i) = ytF(:,end);
 end
 
 % coarse solver
 yc(:,k) = yF(:,k); 
 for i = k+1:n+1
  ytG = Fitz_Time(sG, yc(:,i-1), m); yGc(:,i) = ytG(:,end);
  yc(:,i) = yF(:,i) + yGc(:,i) - yGp(:,i); 
 end
  
  % check errors
 ed = abs(yc(1,end)-yp(1,end));  % relative error
 et = abs(ye(1,end)-yc(1,end));  % true error
 fprintf('iteration %d, %e, %e \n', k-1, ed, et);  

 % check convergence
 if (ed < tol)
  fprintf('sucessful return! \n');
  break
 end

 % update solution
 yGp(:,k:n+1) = yGc(:,k:n+1); 
 yp(:,k:n+1) = yc(:,k:n+1); 
end