In pedroguarderas/CFINI: Computational finance routines

Introduction

1. Loading packages

library( CFINI )
library( plotly )

\begin{equation} \left{ \begin{array}{ll} \dfrac{\partial u}{\partial t} ( t, x ) = \alpha( t, x ) \dfrac{\partial^2 u}{\partial x^2}( t, x ) & \forall ( t, x ) \in [0,T] \times [ y, z ] \ \text{IC:} & \ u( 0, x ) = u_0( x ) & \forall x \in [ y, z ] \ \text{BC:} & \ u( t, y ) = v_1( t ) & \forall t \in [ 0, T ] \ u( t, z ) = v_2( t ) & \forall t \in [ 0, T ] \end{array} \right. \end{equation}

1. Related coefficients

# Diffusion parameter constant
Nt <- 300
Nx <- 150

1. Time grid \begin{equation} t_n = t_{n-1} + n h, \qquad h = \frac{T - 0}{N_t} \end{equation}

t0 <- 0
t1 <- 1
t <- cf_uniform_grid( t0, t1, Nt )

1. Space grid \begin{equation} x_0 = y, \qquad x_i = x_{i-1} + i g,\qquad g = \frac{z - y}{N_x} \end{equation}

y <- -0.5
z <- 0.5
x <- cf_uniform_grid( y, z, Nx )

1. Discretization of parameters \begin{equation} \alpha_{n,i} = \alpha( t_n, x_i ) \end{equation}

alpha <- matrix( 10^(-2.3), Nt, Nx )

1. Initial conditions

uf <- function( x ) if ( abs(x) <= 0.1 ) return( 1 ) else return( 0 )
uf <- Vectorize( uf )
u0 <- sapply( x, FUN = uf )

1. Boundary conditions

v1 <- rep( 0, Nt )
v2 <- rep( 0, Nt )

\begin{equation} A = \begin{bmatrix} b_1 & a_1 & 0 & 0 & \cdots & 0 \ c_1 & b_2 & a_2 & 0 & \cdots & 0 \ 0 & c_2 & b_3 & a_3 & \cdots & 0 \ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \ 0 & 0 & \cdots & c_{n-2} & b_{n-1} & a_{n-1} \ 0 & 0 & \cdots & 0 & c_{n-1} & b_n \end{bmatrix} \end{equation}

1. Euler implicit scheme has the following form for the equation \begin{equation} u_{n+1,i} - u_{n,i} = \alpha_{n,i} \frac{\Delta t_n}{\Delta x_{i}\ \Delta x_{i+1}} ( u_{n+1,i+1} - 2 u_{n+1,i} + u_{n+1,i-1} ) \end{equation}

From previous scheme at every time step $n$ we formulate a tridiagonal problem $A u_n = d$, with the folowing definitions. \begin{eqnarray} \lambda_{n, i} & = & \alpha_{n, i} \frac{\Delta t_n}{\Delta x_{i}\ \Delta x_{i+1}} \ a_i & = & -\lambda_{n,i} \ b_i & = & 1 + 2 \lambda_{n,i} \ c_i & = & -\lambda_{n,i} \ d_i & = & u_{n,i} \end{eqnarray}

Ueu <- cf_diff_solv_euls( alpha, u0, v1, v2, t, x, TRUE )

plot

plot_ly( x = Ueu$t, y = Ueu$x, z = Ueu$u, alpha = 0.8 ) %>% 
  layout( scene = list( xaxis = list(title = "t"),
                        yaxis = list(title = "x"),
                        zaxis = list(title = "u") ) ) %>%
  add_surface()

1. Solving with Crank-Nicolson implicit method \begin{equation} u_{n+1,i} + \theta ( \lambda^2_{n,i} \Delta_x u_{n+1,i+1} - \lambda^1_{n,i} \Delta_x u_{n+1,i} ) = u_{n,i} + ( 1 - \theta ) ( \lambda^2_{n,i} \Delta_x u_{n,i+1} - \lambda^1_{n,i} \Delta_x u_{n,i} ) \end{equation}

\begin{eqnarray} \lambda^1_{n,i} & = & \alpha_{n, i} \frac{\Delta t_n}{\Delta x_{i}\ \Delta x_{i+1}} \ \lambda^2_{n,i} & = & \alpha_{n, i} \frac{\Delta t_n}{\Delta x_{i+1}\ \Delta x_{i+1}} \ a_i & = & -\theta \lambda^1_{n,i} \ b_i & = & 1 + \theta ( \lambda^1_{n,i} + \lambda^2_{n,i} ) \ c_i & = & -\theta \lambda^2_{n,i} \ d_i & = & u_{n,i} + ( 1 - \theta ) ( \lambda^2_{n,i} \Delta_x u_{n,i+1} - \lambda^1_{n,i} \Delta_x u_{n,i} ) \end{eqnarray}

theta <- 0.25
Ucn <- cf_diff_solv_cns( theta, alpha, u0, v1, v2, t, x, TRUE )

plot

plot_ly( x = Ucn$t, y = Ucn$x, z = Ucn$u, alpha = 0.8 ) %>% 
  layout( scene = list( xaxis = list(title = "t"),
                        yaxis = list(title = "x"),
                        zaxis = list(title = "u") ) ) %>%
  add_surface()

1. Plotting results

alph <- max( alpha )
phi <- function( x ) dnorm( x, mean = 0, sd = sqrt( 2 * alph * t[ Nt ] ) )
s <- sapply( x, FUN = function( y ) integrate( function( x, y ) uf( y - x ) * phi( x ), -Inf, Inf, y )$value )

eeu <- matrix( Ueu$u[Nt,] - s, Nx, 1 ) # Euler error
ecn <- matrix( Ucn$u[Nt,] - s, Nx, 1 )

pedroguarderas/CFINI documentation built on Feb. 16, 2024, 2:17 p.m.