From a1d816afe35170ada97cb577989fe59ead2d6dc9 Mon Sep 17 00:00:00 2001 From: Tyler A Date: Wed, 19 May 2021 13:11:12 -0600 Subject: [PATCH] step-78: Black-Scholes Problem. --- doc/doxygen/references.bib | 26 + doc/doxygen/tutorial/tutorial.h.in | 5 + examples/step-78/CMakeLists.txt | 39 ++ examples/step-78/doc/builds-on | 1 + examples/step-78/doc/intro.dox | 363 ++++++++++++ examples/step-78/doc/kind | 1 + examples/step-78/doc/results.dox | 60 ++ examples/step-78/doc/tooltip | 1 + examples/step-78/step-78.cc | 909 +++++++++++++++++++++++++++++ 9 files changed, 1405 insertions(+) create mode 100644 examples/step-78/CMakeLists.txt create mode 100644 examples/step-78/doc/builds-on create mode 100644 examples/step-78/doc/intro.dox create mode 100644 examples/step-78/doc/kind create mode 100644 examples/step-78/doc/results.dox create mode 100644 examples/step-78/doc/tooltip create mode 100644 examples/step-78/step-78.cc diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 1d38ce58b2..fef5c412f1 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1147,6 +1147,32 @@ eprint = {http://dx.doi.org/10.1137/0917003} } +% ------------------------------------ +% Step 78 +% ------------------------------------ + +@article{black1973pricing, + title={The Pricing of Options and Corporate Liabilities}, + author={Black, Fischer and Scholes, Myron}, + journal={The Journal of Political Economy}, + volume={81}, + number={3}, + pages={637--654}, + year={1973} +} + +@article{stoll1969relationship, + title={The relationship between put and call option prices}, + author={Stoll, Hans R}, + journal={The Journal of Finance}, + volume={24}, + number={5}, + pages={801--824}, + year={1969}, + publisher={Wiley Online Library} +} + + % ------------------------------------ % References used elsewhere % ------------------------------------ diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index 4f8b14257a..16bba66822 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -648,6 +648,11 @@ * Interfacing with SUNDIALS' KINSOL nonlinear solver. * * + * + * step-78 + * Solves the Black-Scholes equation for options pricing in 1-D. + * + * * * * diff --git a/examples/step-78/CMakeLists.txt b/examples/step-78/CMakeLists.txt new file mode 100644 index 0000000000..d18ab7b01d --- /dev/null +++ b/examples/step-78/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-78 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-78") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 9.3.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-78/doc/builds-on b/examples/step-78/doc/builds-on new file mode 100644 index 0000000000..1aabbdfdd6 --- /dev/null +++ b/examples/step-78/doc/builds-on @@ -0,0 +1 @@ +step-26 diff --git a/examples/step-78/doc/intro.dox b/examples/step-78/doc/intro.dox new file mode 100644 index 0000000000..087e5c93fd --- /dev/null +++ b/examples/step-78/doc/intro.dox @@ -0,0 +1,363 @@ + +

Introduction

+ +The Black-Scholes equation is a partial differential equation that falls a bit +out of the ordinary scheme. It describes what the fair price of a "European +call" stock option is. Without going into too much detail, a stock "option" is +a contract one can buy from a bank that allows me, but not requires me, to buy +a specific stock at a fixed price $K$ at a fixed future time $T$ in the +future. The question one would then want to answer as a buyer of such an +option is "How much do I think such a contract is worth?", or as the seller +"How much do I need to charge for this contract?", both as a function of the +time $tParticularities of the equation system + +There are a number of oddities in this equation that are worth discussing before +moving on to its numerical solution. First, the "spatial" domain +$\Omega\subset\mathbb{R}$ is unbounded, and thus $S$ can be unbounded in value. +This is because there may be a practical upper bound for stock prices, but not a +conceptual one. The boundary conditions $V(S,t)\rightarrow S$ as +$S\rightarrow \infty$ can then be interpreted as follows: What is the value of +an option that allows me to buy a stock at price $K$ if the stock price (today +or at time $t=T$) is $S\gg K$? One would expect that it is $V\approx S-K$ plus +some adjustment for inflation, or, if we really truly consider huge values of +$S$, we can neglect $K$ and arrive at the statement that the boundary values at +the infinite boundary should be of the form $V\rightarrow S$ as stated above. + +In practice, for us to use a finite element method to solve this, we are going +to need to bound $\Omega$. Since this equation describes prices, and it doesn't +make sense to talk about prices being negative, we will set the lower bound of +$\Omega$ to be 0. Then, for an upper bound, we will choose a very large number, +one that $S$ is not very likely to ever get to. We will call this $S_\text{max}$ +. So, $\Omega=[0,S_\text{max}]$. + +Second, after truncating the domain, we need to ask what boundary values we +should pose at this now finite boundary. To take care of this, we use "put-call" +parity @cite stoll1969relationship. A "pull option" is one in which I am +allowed, but not required, to *sell* a stock at price $K$ to someone at a future +time $T$. This says +@f{align*}{ + V(S,t)+Ke^{-r(T-t)}=P(S,t)+S +@f} +where $V(S,t)$ is the value of the call option, and $P(S,t)$ is the value of the +put option. Since we expect $P(S,t) \rightarrow 0$ as $S \rightarrow \infty$, +this says +@f{align*}{ + V(S,t) \rightarrow S-Ke^{-r(T-t)}, +@f} +and we can use this as a reasonable boundary condition at our finite point +$S_\text{max}$. + +The second complication of the Block-Scholes equation is that we are given a +final condition, and not an initial condition. This is because we know what the +option is worth at time $t=T$: If the stock price at $T$ is $SK$, then I can buy my stock at price $K$ via the option +and immediately sell it again on the market for price $S$, giving me a profit of +$S-K$. In other words, $V(S,T)=S-K$ for $S>K$. So, we only know +values for $V$ at the *end time* but not the initial time -- in fact, finding +out what a fair price at the current time (conventionally taken to be $t=0$) is +what solving these equations is all about. + +This means that this is not an equation that is posed going forward in +time, but in fact going *backward* in time. Thus it makes sense to solve this +problem in reverse by making the change of variables $\tau=T-t$ where now $\tau$ +denotes "time before the strike time $T$". + +With all of this, we finally end up with the following problem: +@f{align*}{ + &-\frac{\partial V}{\partial \tau} + \frac{\sigma^2S^2}{2} \ + \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV=0\ + , \quad\quad &&\forall S\in [0,S_\text{max}], \tau \in [0,T] + \\ + &V(0,\tau) = 0, \ + &&\forall \tau \in [0,T] + \\ + &V(S_\text{max},\tau)=S_\text{max}-Ke^{-r\tau}, \ + &&\forall \tau \in [0,T] + \\ + &V(S,0) = \max(S-K,0) \ + &&\forall S \in [0,S_\text{max}] +@f} + +Conceptually, this is an advection-diffusion-reaction problem for the variable +$V$: There is both a second-order derivative diffusion term, a first-order +derivative advection term, and a zeroth-order reaction term. +We can expect this problem to be a little bit forgiving in practice because for +realistic values of the coefficients, it is diffusive dominated. But, because of +the advective terms in the problem, we will have to be careful with mesh +refinement and time step choice. There is also the issue that the diffusion term + is written in a non-conservative form and so integration by parts is not + immediately obvious. This will be discussed in the next section. + +

Scheme for the numerical solution

+ +We will solve this problem using the fractional step method (of which the +Crank-Nicolson method is a special case with $\theta=\frac 12$; the explicit +Euler method corresponds to $\theta=0$ and the implicit Euler method to +$\theta=1$). So, we first discretize in time, where we would like $V^n(S)$ to +approximate $V(S,\tau_n)$: +@f{align*}{ + 0=&-\frac{V^n(S)-V^{n-1}(S)}{k_n} \\ + &+\frac{\sigma^2S^2}{2}\left[(1-\theta)\frac{d^2V^{n-1}(S)}{dS^2} + \ + \theta \frac{d^2V^{n}(S)}{dS^2}\right] \\ + &+rS\left[(1-\theta)\frac{dV^{n-1}(S)}{dS} + \ + \theta\frac{dV^{n}(S)}{dS}\right] \\ + &-r\left[(1-\theta)V^{n-1}(S) + \theta V^n(S)\right] +@f} +Here, $k_n=\tau_n-\tau_{n-1}$ is the time step size. Given this time +discretization, we can proceed to discretize space by multiplying with test +functions and then integrating by parts. Because there are some interesting +details in this due to the advective and non-advective terms in this equation, +this process will be explained in detail. + +So, we begin by multiplying by test functions, $\{\phi_i(S)\}_{i\in\mathbb{N}}$: +@f{align*}{ + 0=&-\int_0^{S_\text{max}}\phi_i(S)\left[V^n(S)-V^{n-1}(S)\right]dS \\ + &+k_n\int_0^{S_\text{max}}\phi_i(S)\left[\frac{\sigma^2S^2}{2} \ + \left[(1-\theta)\frac{d^2V^{n-1}(S)}{dS^2} + \ + \theta \frac{d^2V^{n}(S)}{dS^2}\right]\right]dS \\ + &+k_n\int_0^{S_\text{max}}\phi_i(S)\left[rS\left[(1-\theta) + \frac{dV^{n-1}(S)}{dS}\ + + \theta\frac{dV^{n}(S)}{dS}\right]\right]dS \\ + &-k_n\int_0^{S_\text{max}}\phi_i(S)\left[r\left[(1-\theta)V^{n-1}(S)\ + + \theta V^n(S)\right]\right]dS +@f} + + +As usual, (1) becomes $-\textbf{M}V^n+\textbf{M}V^{n-1}$ and (4) becomes +$k_n\left[-r(1-\theta)\textbf{M}V^{n-1} - \theta r\textbf{M}V^n\right]$, where +$\textbf{M}_{i,j}=\left(\phi_i(S),\phi_j(S)\right)$, and where we have taken the +liberty of denoting by $V$ not only the function $V(S)$ but also the vector of +nodal values after discretization. + +The interesting parts come from (2) and (3). + + +For (2), we have: +@f{align*}{ + &k_n\int_0^{S_\text{max}}\phi_i(S)\left[\frac{\sigma^2S^2}{2} \ + \left[(1-\theta)\frac{d^2V^{n-1}(S)}{dS^2} + \ + \theta \frac{d^2V^{n}(S)}{dS^2}\right]\right]dS \\ + &=k_n(1-\theta)\int_0^{S_\text{max}}\phi_i(S)\frac{\sigma^2S^2}{2} \ + \frac{d^2V^{n-1}(S)}{dS^2} \ + +k_n\theta\int_0^{S_\text{max}}\phi_i(S)\frac{\sigma^2S^2}{2} \ + \frac{d^2V^{n}(S)}{dS^2} +@f} + +There are two integrals here, that are more or less the same, with the +differences being a slightly different coefficient in front of the integral, +and a different time step for V. Therefore, we will outline this integral in the +general case, and account for the differences at the end. So, consider the +general integral, which we will solve using integration by parts: +@f{align*}{ + &\int_{0}^{S_\text{max}} \phi_i(S)\frac{\sigma^2S^2}{2} + \frac{d^2V^n(S)}{dS^2}dS \\ + &= \phi_i(S)\frac{1}{2}\sigma^2S^2\frac{dV^n(S)}{dS}\Bigg|_0^{S_{max}} - \ + \int_0^{S_\text{max}} \phi_i(S)\sigma^2S\frac{dV^n(S)}{dS}dS - \ + \int_0^{S_\text{max}} \frac{d\phi_i(S)}{dS}\frac{1}{2}\sigma^2S^2 \ + \frac{dV^n(S)}{dS}dS \\ + &= -\int_0^{S_\text{max}} \phi_i(S)\sigma^2S\frac{dV^n(S)}{dS}dS - \ + \int_0^{S_\text{max}} \frac{d\phi_i(S)}{dS}\frac{1}{2}\sigma^2S^2 \ + \frac{dV^n(S)}{dS}dS \\ + &= -\int_0^{S_\text{max}} \phi_i(S)\sigma^2S \sum_j V_j^n + \frac{d\phi_j(S)}{dS}dS \ + -\int_0^{S_\text{max}} \frac{d\phi_i(S)}{dS}\frac{1}{2} \ + \sigma^2S^2 \sum_k V_k^n \frac{d\phi_k(S)}{dS}dS \\ + &= -\sum_j \sigma^2 \int_0^{S_\text{max}} \phi_i(S)S + \frac{d\phi_j(S)}{dS}dS V_j^n\ + - \sum_k \frac{1}{2}\sigma^2 \int_0^{S_\text{max}} \frac{d\phi_i(S)}{dS}S^2\ + \frac{d\phi_k}{dS}dS V_k^n \\ + &= -\sum_j \sigma^2 \left(\phi_i(S)S, \frac{d\phi_j(S)}{dS}\right) V_j^n \ + - \sum_k \frac{1}{2}\sigma^2 \left(\frac{d\phi_i(S)}{dS}S^2,\ + \frac{d\phi_k(S)}{dS}\right) V_k^n \\ + &= -\sigma^2\textbf{B}V^n - \frac{1}{2}\sigma^2\textbf{D}V^n, \quad\quad \ + \textbf{B}_{i,j} = \left(\phi_i(S)S, \frac{d\phi_j(S)}{dS}\right),\ + \textbf{D}_{i,j} = \left(\frac{d\phi_i(S)}{dS}S^2,\frac{d\phi_j(S)}{dS}\right) +@f} + +So, after adding in the constants and exchanging $V^n$ for $V^{n-1}$ where +applicable, we arrive at the following for (2): +@f{align*}{ + &k_n\int_0^{S_\text{max}}\phi_i(S)\left[\frac{\sigma^2S^2}{2} + \left[(1-\theta)\ + \frac{d^2V^{n-1}(S)}{dS^2} + \ + \theta \frac{d^2V^{n}(S)}{dS^2}\right]\right]dS \\ + &= k_n\left[-(1-\theta)\sigma^2\textbf{B}V^{n-1}\ + -(1-\theta)\frac{1}{2}\sigma^2\textbf{D}V^{n-1} \ + -\theta\sigma^2\textbf{B}V^{n} + -\theta\frac{1}{2}\sigma^2\textbf{D}V^{n}\right] +@f} +But, because the matrix $\textbf{B}$ involves an advective term, we will choose +$\theta=0$ there -- in other words, we use an explicit Euler method to treat +advection. Conversely, since the matrix $\textbf{D}$ involves the diffusive term +, we will choose $\theta=1/2$ there -- i.e., we treat diffusion using the second +order Crank-Nicolson method. + +So, we arrive at the following: +@f{align*}{ + k_n\left[-\frac{1}{4}\sigma^2\textbf{D}V^{n-1} \ + -\frac{1}{4}\sigma^2\textbf{D}V^n \ + - \sigma^2\textbf{B}V^{n-1}\right] +@f} + +Now, to handle (3). For this, we will again proceed by considering the general +case like above. + +@f{align*}{ + &\int_{0}^{S_\text{max}} \phi_i(S)rS\frac{dV^n}{dS}dS \\ + &= \phi_i(S)rSV^n\Bigg|_0^{S_\text{max}} - \int_0^{S_\text{max}} + \left[r\phi_i(S) \ + + r\frac{d\phi_i(S)}{dS}S \right]V^ndS \\ + &= -\int_0^{S_\text{max}} r\phi_i(S)V^ndS - \ + \int_0^{S_\text{max}} r\frac{d\phi_i(S)}{dS}SV^ndS \\ + &= -\int_0^{S_\text{max}} r\phi_i(S) \sum_j V_j^n\phi_j(S)dS \ + -\int_0^{S_\text{max}} rS\frac{d\phi_i(S)}{dS} \sum_k V_k\phi_k(S)dS \\ + &= -\sum_j r\left(\phi_i(S), \phi_j(S)\right) V_j^n -\ + \sum_k r\left(S\frac{d\phi_i(S)}{dS}, \phi_k(S)\right)V_k^n \\ + &= -r\textbf{M}V^n -r\textbf{A}V^n, \quad\quad\ + \textbf{M}_{i,j} = \left(\phi_i(S), \phi_j(S)\right),\ + \textbf{A}_{i,j} = \left(S\frac{d\phi_i(S)}{dS}, \phi_j(S)\right) +@f} + +So, again after adding in the constants and exchanging $V^n$ for $V^{n-1}$ where +applicable, we arrive at the following for (3): +@f{align*}{ + &k_n\int_0^{S_\text{max}}\phi_i(S)\left[rS\left[(1-\theta) + \frac{dV^{n-1}(S)}{dS} +\ + \theta\frac{dV^{n}(S)}{dS}\right]\right]dS \\ + &= k_n\left[-(1-\theta)r\textbf{M}V^{n-1} -(1-\theta)r\textbf{A}V^{n-1}\ + -\theta r\textbf{M}V^n -\theta r\textbf{A}V^n\right] +@f} +Just as before, we will use $\theta=0$ for the matrix $\textbf{A}$ and +$\theta=\frac{1}{2}$ for the matrix $\textbf{M}$. So, we arrive at the +following for (3): +@f{align*}{ + k_n\left[-\frac{1}{2}r\textbf{M}V^{n-1} - \frac{1}{2}r\textbf{M}V^n \ + -r\textbf{A}V^{n-1}\right] +@f} + +Now, putting everything together, we obtain the following discrete form for the +Black-Scholes Equation: +@f{align*}{ + 0&= \\ + &-\textbf{M}V^n+\textbf{M}V^{n-1} \\ + & +k_n\left[-\frac{1}{4}\sigma^2\textbf{D}V^{n-1} \ + -\frac{1}{4}\sigma^2\textbf{D}V^n \ + - \sigma^2\textbf{B}V^n \ + -\frac{1}{2}r\textbf{M}V^{n-1} - \frac{1}{2}r\textbf{M}V^n \ + -r\textbf{A}V^n \ + -r\frac{1}{2}\textbf{M}V^{n-1} - \frac{1}{2} r\textbf{M}V^n\right] \\ + &= -\textbf{M}V^n + \textbf{M}V^{n-1} +\ + k_n\left[- \frac{1}{4}\sigma^2\textbf{D}V^{n-1} -\ + \frac{1}{4}\sigma^2\textbf{D}V^n - r\textbf{M}V^{n-1} -\ + r\textbf{M}V^n - \sigma^2\textbf{B}V^{n-1} - r\textbf{A}V^{n-1}\right] +@f} +So, altogether we have: + +@f{equation}{ + 0 = \textbf{M}V^n - \textbf{M}V^{n-1} +\ + k_n\left[ \frac{1}{4}\sigma^2\textbf{D}V^{n-1} +\ + \frac{1}{4}\sigma^2\textbf{D}V^n + r\textbf{M}V^{n-1} + r\textbf{M}V^n +\ + \sigma^2\textbf{B}V^{n-1} + r\textbf{A}V^{n-1}\right]\tag{*} +@f} + +As usual, we can write this with the unknown quantities on the left and the +known ones on the right. This leads to the following linear system that would +have to be solved in each time step: + +@f{align*}{ + \left[\textbf{M}+\frac{1}{4}k_n\sigma^2\textbf{D}+k_nr\textbf{M}\right]V^n\ + =\ + \left[-\frac{1}{4}k_n\sigma^2\textbf{D}-\ + k_nr\textbf{M}+k_n\sigma^2\textbf{B}-\ + k_nr\textbf{A}+\textbf{M}\right]V^{n-1} +@f} + + + + +

Test Case

+For this program, we will use the Method of Manufactured Solutions (MMS) to test + that it is working correctly. This means that we will choose our solution to be + a certain function similar to step-7. For our case, we will use: +@f{align*}{ + V(S,\tau) = -\tau^2 - S^2 + 6\tag{**} +@f} +This means that, using our PDE, we arrive at the following problem: +@f{align*}{ + &-\frac{\partial V}{\partial \tau} +\ + \frac{\sigma^2S^2}{2}\frac{\partial^2 V}{\partial S^2} +\ + rS\frac{\partial V}{\partial S} - rV = f(S,\tau) \\ + &V(0,\tau) = -\tau^2 + 6 \\ + &V(S_\text{max}, \tau) = -S_\text{max}^2 - \tau^2 + 6 \\ + &V(S, 0) = -S^2 + 6 +@f} +Where, $f(S,\tau) = 2\tau - \sigma^2S^2 - 2rS^2 - r(-\tau^2 - S^2 + 6)$. +This set-up now has right hand sides for the equation itself and for the +boundary conditions at $S=0$ that we did not have before, along with "final" +conditions (or, with $\tau$-time "initial conditions") that do not match the +real situation. We will implement this in such a way in the code that it is easy +to exchange -- the introduction of the changes above is just meant to enable the + use of a manufactured solution. + +If the program is working correctly, then it should produce (**) as the +solution. This does mean that we need to modify our variational form somewhat to +account for the non-zero right hand side. + +First, we define the following: +@f{align*}{ + F^n_i = \left(\phi_i(S), f^n(S)\right), && \text{where } f^n(S) =\ + f(S,\tau_n) +@f} +So, we arrive at the new equation: + +@f{align*}{ + \left[\textbf{M}+\frac{1}{4}k_n\sigma^2\textbf{D}+k_nr\textbf{M}\right]V^n\ + =\ + \left[-\frac{1}{4}k_n\sigma^2\textbf{D}-\ + k_nr\textbf{M}+k_n\sigma^2\textbf{B}-\ + k_nr\textbf{A}+\textbf{M}\right]V^{n-1} -\ + k_n\left[\frac{1}{2}F^{n-1}+\frac{1}{2}F^n\right] +@f} + +We then solve this equation as outlined above. diff --git a/examples/step-78/doc/kind b/examples/step-78/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-78/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-78/doc/results.dox b/examples/step-78/doc/results.dox new file mode 100644 index 0000000000..db4bb2d2ad --- /dev/null +++ b/examples/step-78/doc/results.dox @@ -0,0 +1,60 @@ +

Results

+ + +Below is the output of the program: +@code +=========================================== +Number of active cells: 1 +Number of degrees of freedom: 2 + +Time step 1 at t=0.0002 + 0 CG iterations. +[...] +=========================================== +Number of active cells: 128 +Number of degrees of freedom: 129 + +Time step 1 at t=0.0002 + 2 CG iterations. +[...] +Time step 5001 at t=1.0002 + 5 CG iterations. +Cycle 7: + Number of active cells: 128 + Number of degrees of freedom: 129 + +cells dofs L2 H1 Linfty + 1 2 1.667e-01 5.774e-01 2.222e-01 + 2 3 3.906e-02 2.889e-01 5.380e-02 + 4 5 9.679e-03 1.444e-01 1.357e-02 + 8 9 2.405e-03 7.218e-02 3.419e-03 + 16 17 5.967e-04 3.609e-02 8.597e-04 + 32 33 1.457e-04 1.804e-02 2.155e-04 + 64 65 3.306e-05 9.022e-03 5.388e-05 + 128 129 5.014e-06 4.511e-03 1.342e-05 + +n cells H1 L2 + 1 5.774e-01 - - 1.667e-01 - - + 2 2.889e-01 2.00 1.00 3.906e-02 4.27 2.09 + 4 1.444e-01 2.00 1.00 9.679e-03 4.04 2.01 + 8 7.218e-02 2.00 1.00 2.405e-03 4.02 2.01 + 16 3.609e-02 2.00 1.00 5.967e-04 4.03 2.01 + 32 1.804e-02 2.00 1.00 1.457e-04 4.10 2.03 + 64 9.022e-03 2.00 1.00 3.306e-05 4.41 2.14 + 128 4.511e-03 2.00 1.00 5.014e-06 6.59 2.72 + +@endcode + +What is more interesting is the output of the convergence tables. They are +outputted into the console, as well into a LaTex file. The convergence tables +are shown above. Here, you can see that the the solution has a convergence rate +of $\mathcal{O}(h)$ with respect to the $H^1$-norm, and the solution has a convergence rate +of $\mathcal{O}(h^2)$ with respect to the $L^2$-norm. + + +Below is the visualization of the solution. + +
+ Solution of the MMS problem. +
diff --git a/examples/step-78/doc/tooltip b/examples/step-78/doc/tooltip new file mode 100644 index 0000000000..c16d4734c6 --- /dev/null +++ b/examples/step-78/doc/tooltip @@ -0,0 +1 @@ +Black-Scholes equation for stock options. diff --git a/examples/step-78/step-78.cc b/examples/step-78/step-78.cc new file mode 100644 index 0000000000..431e50cbe5 --- /dev/null +++ b/examples/step-78/step-78.cc @@ -0,0 +1,909 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2021 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Author: Tyler Anderson, Colorado State University, 2021 + */ + + +// @sect3{Include files} + +// The program starts with the usual include files, all of which you should have +// seen before by now: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +// Then the usual placing of all content of this program into a namespace and +// the importation of the deal.II namespace into the one we will work in. We +// also define an identifier to allow for the MMS code to be run when +// MMS is defined. Otherwise, the program solves the original +// problem: +namespace BlackScholesSolver +{ + using namespace dealii; + +#define MMS + + // @sect3{Solution Class} + + // This section creates a class for the known solution when testing using the + // MMS. Here I am using $v(\tau,S) = -\tau^2 -S^2 + 6$ for my solution. We + // need to include the solution equation and the gradient for the H1 seminorm + // calculation. + template + class Solution : public Function + { + public: + Solution(const double maturity_time); + + virtual double value(const Point & p, + const unsigned int component = 0) const override; + + virtual Tensor<1, dim> + gradient(const Point & p, + const unsigned int component = 0) const override; + + private: + const double maturity_time; + }; + + + template + Solution::Solution(const double maturity_time) + : maturity_time(maturity_time) + { + Assert(dim == 1, ExcNotImplemented()); + } + + + template + double Solution::value(const Point & p, + const unsigned int component) const + { + return -Utilities::fixed_power<2, double>(p(component)) - + Utilities::fixed_power<2, double>(this->get_time()) + 6; + } + + + template + Tensor<1, dim> Solution::gradient(const Point & p, + const unsigned int component) const + { + return Point(-2 * p(component)); + } + + + + // @sect3{Equation Data} + + // In the following classes and functions, we implement the right hand side + // and boundary values that define this problem and for which we need function + // objects. The right hand side is chosen as discussed at the end of the + // introduction. + // + // First, we handle the initial condition. + template + class InitialConditions : public Function + { + public: + InitialConditions(const double strike_price); + + virtual double value(const Point & p, + const unsigned int component = 0) const override; + + private: + const double strike_price; + }; + + + template + InitialConditions::InitialConditions(const double strike_price) + : strike_price(strike_price) + {} + + + template + double InitialConditions::value(const Point & p, + const unsigned int component) const + { +#ifdef MMS + return -Utilities::fixed_power<2, double>(p(component)) + 6; +#else + return std::max(p(component) - strike_price, 0.); +#endif + } + + + + // Next, we handle the left boundary condition. + template + class LeftBoundaryValues : public Function + { + public: + virtual double value(const Point & p, + const unsigned int component = 0) const override; + }; + + + template + double LeftBoundaryValues::value(const Point &, + const unsigned int /*component*/) const + { +#ifdef MMS + return -Utilities::fixed_power<2, double>(this->get_time()) + 6; +#else + return 0.; +#endif + } + + + + // Then, we handle the right boundary condition. + template + class RightBoundaryValues : public Function + { + public: + RightBoundaryValues(const double strike_price, const double interest_rate); + + virtual double value(const Point & p, + const unsigned int component = 0) const override; + + private: + const double strike_price; + const double interest_rate; + }; + + + template + RightBoundaryValues::RightBoundaryValues(const double strike_price, + const double interest_rate) + : strike_price(strike_price) + , interest_rate(interest_rate) + {} + + + template + double RightBoundaryValues::value(const Point & p, + const unsigned int component) const + { +#ifdef MMS + return -Utilities::fixed_power<2, double>(p(component)) - + Utilities::fixed_power<2, double>(this->get_time()) + 6; +#else + return (p(component) - strike_price) * + exp((-interest_rate) * (this->get_time())); +#endif + } + + + + // Finally, we handle the right hand side. + template + class RightHandSide : public Function + { + public: + RightHandSide(const double asset_volatility, const double interest_rate); + + virtual double value(const Point & p, + const unsigned int component = 0) const override; + + private: + const double asset_volatility; + const double interest_rate; + }; + + + template + RightHandSide::RightHandSide(const double asset_volatility, + const double interest_rate) + : asset_volatility(asset_volatility) + , interest_rate(interest_rate) + {} + + + template + double RightHandSide::value(const Point & p, + const unsigned int component) const + { +#ifdef MMS + return 2 * (this->get_time()) - + Utilities::fixed_power<2, double>(asset_volatility * p(component)) - + 2 * interest_rate * Utilities::fixed_power<2, double>(p(component)) - + interest_rate * + (-Utilities::fixed_power<2, double>(p(component)) - + Utilities::fixed_power<2, double>(this->get_time()) + 6); +#else + (void)p; + (void)component; + return 0.0; +#endif + } + + + + // @sect3{The BlackScholes Class} + + // The next piece is the declaration of the main class of this program. This + // is very similar to the Step-26 tutorial, with some modifications. New + // matrices had to be added to calculate the A and B matrices, as well as the + // $V_{diff}$ vector mentioned in the introduction. We also define the + // parameters used in the problem. + // + // - maximum_stock_price: The imposed upper bound on the spatial + // domain. This is the maximum allowed stock price. + // - maturity_time: The upper bound on the time domain. This is + // when the option expires.\n + // - asset_volatility: The volatility of the stock price.\n + // - interest_rate: The risk free interest rate.\n + // - strike_price: The aggreed upon price that the buyer will + // have the option of purchasing the stocks at the expiration time. + // + // Some slight differences between this program and step-26 are the creation + // of the a_matrix and the b_matrix, which is + // described in the introduction. We then also need to store the current time, + // the size of the time step, and the number of the current time step. + // Next, we will store the output into a DataOutStack + // variable because we will be layering the solution at each time on top of + // one another to create the solution manifold. Then, we have a variable that + // stores the current cycle and number of cycles that we will run when + // calculating the solution. The cycle is one full solution calculation given + // a mesh. We refine the mesh once in between each cycle to exhibit the + // convergence properties of our program. Finally, we store the convergence + // data into a convergence table. + // + // As far as member functions are concerned, we have a function that + // calculates the convergence information for each cycle, called + // process_solution. This is just like what is done in step-7. + template + class BlackScholes + { + public: + BlackScholes(); + + void run(); + + private: + void setup_system(); + void solve_time_step(); + void refine_grid(); + void process_solution(); + void add_results_for_output(); + void write_convergence_table(); + + const double maximum_stock_price; + const double maturity_time; + const double asset_volatility; + const double interest_rate; + const double strike_price; + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix mass_matrix; + SparseMatrix laplace_matrix; + SparseMatrix a_matrix; + SparseMatrix b_matrix; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + double time; + double time_step; + unsigned int timestep_number; + + const double theta; + const unsigned int n_cycles; + + DataOutStack data_out_stack; + std::vector solution_names; + + ConvergenceTable convergence_table; + }; + + // @sect3{The BlackScholes Implementation} + + // Now, we get to the implementation of the main class. We will set the values + // for the various parameters used in the problem. These were chosen because + // they are fairly normal values for these parameters. Although the stock + // price has no upper bound in reality (it is in fact infinite), we impose + // an upper bound that is twice the strike price. This is a somewhat arbitrary + // choice to be twice the strike price, but it is large enought to see the + // interesting parts of the solution. + template + BlackScholes::BlackScholes() + : maximum_stock_price(1.) + , maturity_time(1.) + , asset_volatility(.2) + , interest_rate(0.05) + , strike_price(0.5) + , fe(1) + , dof_handler(triangulation) + , time(0.0) + , timestep_number(0) + , theta(0.5) + , n_cycles(3) + { + Assert(dim == 1, ExcNotImplemented()); + } + + // @sect4{BlackScholes::setup_system} + + // The next function sets up the DoFHandler object, computes + // the constraints, and sets the linear algebra objects to their correct + // sizes. We also compute the mass matrix here by calling a function from the + // library. We will compute the other 3 matrices next, because these need to + // be computed 'by hand'. + // + // Note, that the time step is initialized here because the maturity time was + // needed to compute the time step. + template + void BlackScholes::setup_system() + { + dof_handler.distribute_dofs(fe); + + time_step = maturity_time / 5000.; + + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + constraints, + /*keep_constrained_dofs = */ true); + sparsity_pattern.copy_from(dsp); + + mass_matrix.reinit(sparsity_pattern); + laplace_matrix.reinit(sparsity_pattern); + a_matrix.reinit(sparsity_pattern); + b_matrix.reinit(sparsity_pattern); + system_matrix.reinit(sparsity_pattern); + + MatrixCreator::create_mass_matrix(dof_handler, + QGauss(fe.degree + 1), + mass_matrix); + + // Below is the code to create the Laplace matrix with non-constant + // coefficients. This corresponds to the matrix D in the introduction. This + // non-constant coefficient is represented in the + // current_coefficient variable. + const unsigned int dofs_per_cell = fe.dofs_per_cell; + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + QGauss quadrature_formula(fe.degree + 1); + FEValues fe_values(fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + std::vector local_dof_indices(dofs_per_cell); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell_matrix = 0.; + fe_values.reinit(cell); + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + const double current_coefficient = + fe_values.quadrature_point(q_index).square(); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + cell_matrix(i, j) += + (current_coefficient * // (x_q)^2 + fe_values.shape_grad(i, q_index) * // grad phi_i(x_q) + fe_values.shape_grad(j, q_index) * // grad phi_j(x_q) + fe_values.JxW(q_index)); // dx + } + } + cell->get_dof_indices(local_dof_indices); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + laplace_matrix.add(local_dof_indices[i], + local_dof_indices[j], + cell_matrix(i, j)); + } + } + + // Now we will create the A matrix. Below is the code to create the matrix A + // as discussed in the introduction. The non constant coefficient is again + // represented in the current_coefficient variable. + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell_matrix = 0.; + fe_values.reinit(cell); + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + const Tensor<1, dim> current_coefficient = + fe_values.quadrature_point(q_index); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + { + cell_matrix(i, j) += + (current_coefficient * // x_q + fe_values.shape_grad(i, q_index) * // grad phi_i(x_q) + fe_values.shape_value(j, q_index) * // phi_j(x_q) + fe_values.JxW(q_index)); // dx + } + } + } + cell->get_dof_indices(local_dof_indices); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + a_matrix.add(local_dof_indices[i], + local_dof_indices[j], + cell_matrix(i, j)); + } + } + + // Finally we will create the matrix B. Below is the code to create the + // matrix B as discussed in the introduction. The non constant coefficient + // is again represented in the current_coefficient variable. + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell_matrix = 0.; + fe_values.reinit(cell); + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + const Tensor<1, dim> current_coefficient = + fe_values.quadrature_point(q_index); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + cell_matrix(i, j) += + (current_coefficient * // x_q + fe_values.shape_value(i, q_index) * // phi_i(x_q) + fe_values.shape_grad(j, q_index) * // grad phi_j(x_q) + fe_values.JxW(q_index)); // dx + } + } + cell->get_dof_indices(local_dof_indices); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + b_matrix.add(local_dof_indices[i], + local_dof_indices[j], + cell_matrix(i, j)); + } + } + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + // @sect4{BlackScholes::solve_time_step} + + // The next function is the one that solves the actual linear system for a + // single time step. The only interesting thing here is that the matrices + // we have built are symmetric positive definite, so we can use the + // conjugate gradient method. + template + void BlackScholes::solve_time_step() + { + SolverControl solver_control(1000, 1e-12); + SolverCG> cg(solver_control); + PreconditionSSOR> preconditioner; + preconditioner.initialize(system_matrix, 1.0); + cg.solve(system_matrix, solution, system_rhs, preconditioner); + constraints.distribute(solution); + std::cout << " " << solver_control.last_step() << " CG iterations." + << std::endl; + } + + // @sect4{BlackScholes::add_results_for_output} + + // This is simply the function to stitch the solution peices together. For + // this, we create a new layer at each time, and then add the solution vector + // for that timestep. The function then stitches this together with the old + // solutions using 'build_patches'. + template + void BlackScholes::add_results_for_output() + { + data_out_stack.new_parameter_value(time, time_step); + data_out_stack.attach_dof_handler(dof_handler); + data_out_stack.add_data_vector(solution, solution_names); + data_out_stack.build_patches(2); + data_out_stack.finish_parameter_value(); + } + + // @sect4{BlackScholes::refine_grid} + + // It is somewhat unnecessary to have a function for the global refinement + // that we do. The reason for the function is to allow for the possibility of + // an adaptive refinement later. + template + void BlackScholes::refine_grid() + { + triangulation.refine_global(1); + } + + // @sect4{BlackScholes::process_solution} + + // This is where we calculate the convergence and error data to evaluate the + // effectiveness of the program. Here, we calculate the $L^2$, $H^1$ and + // $L^{\infty}$ norms. + template + void BlackScholes::process_solution() + { + Solution sol(maturity_time); + sol.set_time(time); + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(dof_handler, + solution, + sol, + difference_per_cell, + QGauss(fe.degree + 1), + VectorTools::L2_norm); + const double L2_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + VectorTools::integrate_difference(dof_handler, + solution, + sol, + difference_per_cell, + QGauss(fe.degree + 1), + VectorTools::H1_seminorm); + const double H1_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::H1_seminorm); + const QTrapezoid<1> q_trapezoid; + const QIterated q_iterated(q_trapezoid, fe.degree * 2 + 1); + VectorTools::integrate_difference(dof_handler, + solution, + sol, + difference_per_cell, + q_iterated, + VectorTools::Linfty_norm); + const double Linfty_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::Linfty_norm); + const unsigned int n_active_cells = triangulation.n_active_cells(); + const unsigned int n_dofs = dof_handler.n_dofs(); + convergence_table.add_value("cells", n_active_cells); + convergence_table.add_value("dofs", n_dofs); + convergence_table.add_value("L2", L2_error); + convergence_table.add_value("H1", H1_error); + convergence_table.add_value("Linfty", Linfty_error); + } + + //@sect4{BlackScholes::write_convergence_table } + + // This next part is building the convergence and error tables. By this, we + // need to set the settings for how to output the data that was calculated + // during BlackScholes::process_solution. First, we will create + // the headings and set up the cells properly. During this, we will also + // prescribe the precision of our results. Then we will write the calculated + // errors based on the $L^2$, $H^1$, and $L^{\infty}$ norms to the console and + // to the error LaTeX file. + template + void BlackScholes::write_convergence_table() + { + convergence_table.set_precision("L2", 3); + convergence_table.set_precision("H1", 3); + convergence_table.set_precision("Linfty", 3); + convergence_table.set_scientific("L2", true); + convergence_table.set_scientific("H1", true); + convergence_table.set_scientific("Linfty", true); + convergence_table.set_tex_caption("cells", "\\# cells"); + convergence_table.set_tex_caption("dofs", "\\# dofs"); + convergence_table.set_tex_caption("L2", "@f$L^2@f$-error"); + convergence_table.set_tex_caption("H1", "@f$H^1@f$-error"); + convergence_table.set_tex_caption("Linfty", "@f$L^\\infty@f$-error"); + convergence_table.set_tex_format("cells", "r"); + convergence_table.set_tex_format("dofs", "r"); + std::cout << std::endl; + convergence_table.write_text(std::cout); + std::string error_filename = "error"; + error_filename += "-global"; + error_filename += ".tex"; + std::ofstream error_table_file(error_filename); + convergence_table.write_tex(error_table_file); + + // Next, we will make the convergence table. We will again write this to + // the console and to the convergence LaTex file. + convergence_table.add_column_to_supercolumn("cells", "n cells"); + std::vector new_order; + new_order.emplace_back("n cells"); + new_order.emplace_back("H1"); + new_order.emplace_back("L2"); + convergence_table.set_column_order(new_order); + convergence_table.evaluate_convergence_rates( + "L2", ConvergenceTable::reduction_rate); + convergence_table.evaluate_convergence_rates( + "L2", ConvergenceTable::reduction_rate_log2); + convergence_table.evaluate_convergence_rates( + "H1", ConvergenceTable::reduction_rate); + convergence_table.evaluate_convergence_rates( + "H1", ConvergenceTable::reduction_rate_log2); + std::cout << std::endl; + convergence_table.write_text(std::cout); + std::string conv_filename = "convergence"; + conv_filename += "-global"; + switch (fe.degree) + { + case 1: + conv_filename += "-q1"; + break; + case 2: + conv_filename += "-q2"; + break; + default: + Assert(false, ExcNotImplemented()); + } + conv_filename += ".tex"; + std::ofstream table_file(conv_filename); + convergence_table.write_tex(table_file); + } + + // @sect4{BlackScholes::run} + + // Now we get to the main driver of the program. This is where we do all the + // work of looping through the time steps and calculating the solution vector + // each time. Here at the top, we set the initial refinement value and then + // create a mesh. Then we refine this mesh once. Next, we set up the + // data_out_stack object to store our solution. Finally, we start a for loop + // to loop through the cycles. This lets us recalculate a solution for each + // successive mesh refinement. At the beginning of each iteration, we need to + // reset the time and time step number. We introduce an if statement to + // accomplish this because we don't want to do this on the first iteration. + template + void BlackScholes::run() + { + GridGenerator::hyper_cube(triangulation, 0.0, maximum_stock_price, true); + triangulation.refine_global(0); + + solution_names.emplace_back("u"); + data_out_stack.declare_data_vector(solution_names, + DataOutStack::dof_vector); + + Vector vmult_result; + Vector forcing_terms; + + for (unsigned int cycle = 0; cycle < n_cycles; cycle++) + { + if (cycle != 0) + { + refine_grid(); + time = 0.0; + timestep_number = 0; + } + + setup_system(); + + std::cout << std::endl + << "===========================================" << std::endl + << "Cycle " << cycle << ':' << std::endl + << "Number of active cells: " + << triangulation.n_active_cells() << std::endl + << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << std::endl; + + VectorTools::interpolate(dof_handler, + InitialConditions(strike_price), + solution); + + if (cycle == (n_cycles - 1)) + { + add_results_for_output(); + } + + // Next, we run the main loop which runs until we exceed the maturity + // time. We first compute the right hand side of the equation, which is + // described in the introduction. Recall that it contains the term + // $\left[-\frac{1}{4}k_n\sigma^2\mathbf{D}-k_nr\mathbf{M}+k_n\sigma^2 + // \mathbf{B}-k_nr\mathbf{A}+\mathbf{M}\right]V^{n-1}$. We put these + // terms into the variable system_rhs, with the help of a temporary + // vector: + vmult_result.reinit(dof_handler.n_dofs()); + forcing_terms.reinit(dof_handler.n_dofs()); + while (time < maturity_time) + { + time += time_step; + ++timestep_number; + std::cout << "Time step " << timestep_number << " at t=" << time + << std::endl; + + mass_matrix.vmult(system_rhs, solution); + + laplace_matrix.vmult(vmult_result, solution); + system_rhs.add( + (-1) * (1 - theta) * time_step * + Utilities::fixed_power<2, double>(asset_volatility) * 0.5, + vmult_result); + mass_matrix.vmult(vmult_result, solution); + + system_rhs.add((-1) * (1 - theta) * time_step * interest_rate * 2, + vmult_result); + + a_matrix.vmult(vmult_result, solution); + system_rhs.add((-1) * time_step * interest_rate, vmult_result); + + b_matrix.vmult(vmult_result, solution); + system_rhs.add( + (-1) * Utilities::fixed_power<2, double>(asset_volatility) * + time_step * 1, + vmult_result); + + // The second piece is to compute the contributions of the source + // terms. This corresponds to the term $-k_n\left[\frac{1}{2}F^{n-1} + // +\frac{1}{2}F^n\right]$. The following code calls + // VectorTools::create_right_hand_side to compute the vectors $F$, + // where we set the time of the right hand side (source) function + // before we evaluate it. The result of this all ends up in the + // forcing_terms variable: + RightHandSide rhs_function(asset_volatility, interest_rate); + rhs_function.set_time(time); + VectorTools::create_right_hand_side(dof_handler, + QGauss(fe.degree + 1), + rhs_function, + forcing_terms); + forcing_terms *= time_step * theta; + system_rhs -= forcing_terms; + + rhs_function.set_time(time - time_step); + VectorTools::create_right_hand_side(dof_handler, + QGauss(fe.degree + 1), + rhs_function, + forcing_terms); + forcing_terms *= time_step * (1 - theta); + system_rhs -= forcing_terms; + + // Next, we add the forcing terms to the ones that come from the + // time stepping, and also build the matrix $\left[\mathbf{M}+ + // \frac{1}{4}k_n\sigma^2\mathbf{D}+k_nr\mathbf{M}\right]$ that we + // have to invert in each time step. The final piece of these + // operations is to eliminate hanging node constrained degrees of + // freedom from the linear system: + system_matrix.copy_from(mass_matrix); + system_matrix.add( + (theta)*time_step * + Utilities::fixed_power<2, double>(asset_volatility) * 0.5, + laplace_matrix); + system_matrix.add((time_step)*interest_rate * theta * (1 + 1), + mass_matrix); + + constraints.condense(system_matrix, system_rhs); + + // There is one more operation we need to do before we can solve it: + // boundary values. To this end, we create a boundary value object, + // set the proper time to the one of the current time step, and + // evaluate it as we have done many times before. The result is used + // to also set the correct boundary values in the linear system: + { + RightBoundaryValues right_boundary_function(strike_price, + interest_rate); + LeftBoundaryValues left_boundary_function; + right_boundary_function.set_time(time); + left_boundary_function.set_time(time); + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + left_boundary_function, + boundary_values); + VectorTools::interpolate_boundary_values(dof_handler, + 1, + right_boundary_function, + boundary_values); + MatrixTools::apply_boundary_values(boundary_values, + system_matrix, + solution, + system_rhs); + } + + // With this out of the way, all we have to do is solve the system, + // generate graphical data on the last cycle, and create the + // convergence table data. + solve_time_step(); + + if (cycle == (n_cycles - 1)) + { + add_results_for_output(); + } + } +#ifdef MMS + process_solution(); +#endif + } + + const std::string filename = "solution.vtk"; + std::ofstream output(filename); + data_out_stack.write_vtk(output); + +#ifdef MMS + write_convergence_table(); +#endif + } + +} // namespace BlackScholesSolver + +// @sect3{The main Function} + +// Having made it this far, there is, again, nothing much to discuss for the +// main function of this program: it looks like all such functions since step-6. +int main() +{ + try + { + using namespace BlackScholesSolver; + + BlackScholes<1> black_scholes_solver; + black_scholes_solver.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; +} -- 2.39.5