## 2 Results
### Computational Experiments  
Consider (9)-(14). Similar to the discussions by [PERSON] et al. [14], we fix \(X=8,\ Y=1\). We first concentrate on experiments with \(\rho=-0.5\) and \(T=0.5\). Next, to test against extreme cases in the option market, we proceed with \(\rho=-1\) and \(T=5\). For demonstrating the numerical solution and its rate of convergence estimates, we first consider uniform spatial grids. To this end, we denote  
\[h_{m}=h,\ k_{n}=k=\sigma h,\ \ m=1,2,\ldots,M;\ n=1,2,\ldots,N.\]  
Results over nonuniform grids will be presented later. Some key parameters used are shown in Table 1. Further, a Crank-Nicolson type temporal integrator will be utilized for advancing our semi-discretized system (25), (28), with \(\Delta\tau\) as the temporal step [19]. It is known that \(\lambda=\Delta\tau/c^{2}\), where \(c=\min\left\{h,k\right\}\), plays the effective role of the Courant number [21; 22]. We experiment with different values of \(\lambda\), varying from \(0.5\) to \(1\). Our semi-discretized scheme is expected to be up to the first order in convergence in space. To numerically examine this by experiment, we employ a generalized Milne's device [10; 18]. Then, for a selected terminal time \(T\), we denote the numerical solution at point \((x_{m},y_{n},T)\), \(1\leq m\leq M,1\leq n\leq N\), as \(u_{m,\eta,h}\) for any particular spatial step \(0<h\ll 1\). Likewise, we let \begin{table}
\begin{tabular}{c c} \hline \hline
**Key Parameter** & **Value Used** \\ \hline Strike price & \(K=100\) \\ Volatility of volatility & \(\sigma=1\) \\ Risk-free interest rate & \(r=0.05\) \\ Mean reversion speed & \(\kappa=2\) \\ Long-run mean of volatility & \(\eta=0.1\) \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Key parameter values for numerical simulations. \(u_{m,n;k/2}\) and \(u_{m,n;k/4}\) be the computed solutions obtained by using \(h/2\) and \(h/4\), respectively. Thus, the point-wise rate of spatial convergence at \(T\) can be evaluated via  
\[R^{h}_{m,n}\approx\frac{1}{\text{ln}2}\text{ln}\frac{\left|u_{m,n;k}-u_{m,n;k/2 }\right|}{\left|u_{m,n;k/2}-u_{m,n;k/4}\right|}. \tag{32}\]  
Most of our experiments are accomplished on Apple workstations. MATLAB platforms without parallelizations are used throughout our operations. Let \(h=0.01\) and \(\sigma=1\). For simplicity of notation, we use the same letter \(v\) for the approximate solution to (4). We show the solution \(v\) for \(\rho=-0.5\) and \(\rho=-1\) in Figures 3 and 4, respectively. To see more precise solution profiles, we show corresponding contour maps next to the surfaces. It can be observed that the European put option price is a decreasing function of the stock price \(S\). This coincides well with the financial theory that a put option price should have a negative correlation with the underline stock price [11; 23]. To examine further the delicate relationship between a put option price and its volatility, we plot an average numerical solution \(\mathcal{G}(y,t)\) taken across different stock prices with respect to the volatility in Figure 5. The simulated computational result is exactly what we would expect, since a put option price should be positively correlated with volatility [22; 23]. We plot the computed rate of convergence surfaces for the cases \(\rho=-0.5\) and \(\rho=-1\) in Figures 6 and 7, respectively. In addition, a summary of point-wise convergence rates for the circumstance as \(\rho=-0.5,\ T=0.5\) on different spatial grids is given in Table 2. Minor disturbances can be observed in regions where the solution changes fast, particularly in the extreme situations with \(\rho=-1\), as demonstrated in Figure 7. These results are consistent with those from well-established high-order schemes [2; 4; 8; 13; 14]. However, due to the low-order nature of our scheme, to get the same accuracy we have to employ relatively small mesh sizes. This results in longer running  
Figure 4: **(Left)** Price of an European put option at \(T=5\) and for \(\rho=-1\); **(Right)** Corresponding contour map.