有限体积交替方向隐式方法

因此，必须将空间域截断为[xmin，xmax]，其中边界选择得离x足够远，因此运行误差可以忽略不计。回想一下，σ，u的形式可以自然地将PDE的空间域限制为例如x≥ 在后一种情况下，下边界自然定义为xmin=0。如前所述，定义空间网格xmin=x<x<…<xm=xmax，letxi=xi- xi-1网格宽度，带x=xm+1=0，定义-0.5=xi-xi=xi-1+XI代表2≤ 我≤ m、 x0.5=X和xm+0.5=xm。这将生成一个以顶点为中心的网格，其中包含单元格Ohmi=【xi】-0.5，xi+0.5]。设Pi（τ）表示精确单元平均值的近似值spi（τ）=xi+0.5- xi-0.5ZOhmip（x，τ）dx，设P为包含这些近似值的向量。与上一节类似（见方程式（2.4）和（2.5）以及[20]），我们定义了公式pi（τ）=[fi-0.5（P，τ）-fi+0.5（P，τ）]xi+xi+1（2.9），其中数值flux由fi±0.5（P，τ）=fa，i±0.5（P，τ）+fd，i±0.5（P，τ），其中fa，i±0.5（P，τ）≈ u（xi±0.5，τ）p（xi±0.5，τ），（2.10）和fd，i±0.5（p，τ）≈ -x个σ（x，τ）p（x，τ）|x=xi±0.5。（2.11）为了便于表述，从现在起，我们省略参数对τ的依赖性，并将ui±0.5=u（xi±0.5，τ）和σi=σ（xi，τ）。注意，f0,5（P，τ），分别为fm+0.5（P，τ），对应于边界xmin=x处的flux，分别为xmax=xm。PDE（2.8）的平流部分以保守形式书写。对于内部单元边界，即xi-0.5带2≤ 我≤ m、 我们考虑二阶中央FV方案，参见[20]，并定义fa，i-0.5（P，τ）in（2.10）asfa，i-0.5（P，τ）=ui-0.5Pi-1（τ）+π（τ）。差异部分并非以保守形式编写，因此不可能将标准FV方案直接应用于该术语。（2.11）的二阶FV方案的想法，参见示例。

[20] ，由definingfd，i概括-0.5（P，τ）=-σiPi（τ）-σi-1Pi-1（τ）xifor 2≤ 我≤ m、 很容易看出σip（xi，τ）-σi-1p（xi-1，τ）xi是的二阶近似值xi点的x（σp）-0.5这解释了这种离散化的选择。将这些表达式插回（2.9）中，我们得到pi（τ）=σi-1Pi-1（τ）xi(xi+xi+1）-σiPi（τ）xixi+1+σi+1Pi+1（τ）xi+1(xi+xi+1）（2.12）+ui-0.5Pi-1（τ）+π（τ）- ui+0.5Pi（τ）+Pi+1（τ）xi+xi+1，用于2≤ 我≤ m级- 请注意，通过在非均匀空间网格上应用二阶中心FD格式，参见[14]，在x（σp），一个将以相同的离散化结束。为了完成这种半离散化，还必须在截断域的边界处进行定义，以确保总质量守恒。假设p代表密度函数，那么z∞-∞p（x，τ）dx=1，τ>0，henceZ∞-∞τpdx=Z∞-∞h类x个σp-x（up）idx=0。假设[xmin，xmax]选择得足够宽，上述条件可近似为x个σp- （up）x=xmaxx=xmin=0，反映了间隔[xmin，xmax]上的总流量为零的事实。如上所述，对于某些系数函数的选择，空间域自然受到限制。例如，如果thePDE（2.8）源于非负过程，则xminca应设置为零。上述无flux边界条件仍然适用于自然限制域。从理论上讲，可能在其中一个边界处存在正flux，而在另一个边界处存在完全相同的负flux。然而，由于该解表示密度函数，因此更现实的做法是，在每个边界处都没有质量进出。

有鉴于此，我们假设以下边界条件成立：x个σp- （up）x=xmin=0，（2.13）x个σp- （up）x=xmax=0。第一个条件的数值等价物是，x=x最小值时的flux为零，即f（P，τ）≡ f0.5（P，τ）=0。这可以通过创建重影点x=x来实现-X并使用（2.13）通过σP（τ）确定重影点处σP（τ）的值-σP（τ）2x个- uP（τ）=0，其中我们利用了单元平均值形成点值的二阶近似值这一事实。转向（2.11），我们定义了不同的流量fd，0.5at xasfd，0.5（P，τ）=-σP（τ）-σP（τ）2x=-uP（τ）。由于xis是第一个单元格的左边界，因此从前倾部分开始的边界x上的flux（见（2.10））可以近似为fa，0.5=uP（τ）。将这些表达式插入（2.9）中，我们得到p（τ）=-f1.5（P，τ）x+x=-f1.5（P，τ）x、 （2.14）可以类似地处理xmax处的边界条件，以获得pm（τ）=fm-0.5（P，τ）xm。（2.15）通过以这种方式执行边界条件的离散化，可以得出f0.5（P，τ）=fm+0.5（P，τ），并且我们确保质量在数值格式中守恒。结合（2.12）、（2.14）和（2.15），我们可以看到，对于给定矩阵A（τ），对于τ>0，总离散化可以写成ODEP（τ）=A（τ）P（τ）（2.16）的系统。由于Pi（τ）代表细胞平均值，因此很自然地定义初始向量asPi（0）=xi+xi+1if X∈ 【xi】-0.5，xi+0.5]，否则为0。一般来说，常微分方程组（2.16）的精确解不能用解析方法计算，只能依靠数值方法来近似。由于上述离散化往往导致Stiff半离散系统，因此广泛考虑了合适的隐式时间步进格式，如Crank–Nicolson格式，参见示例【29】。2.3数值试验在本小节中，通过考虑两个实际示例来测试FV离散化的性能。

作为第一个示例，考虑SDEdSτ=（rd- rf）Sτdτ+σBSSτdWτ，σBS>0，对应于经典的Black-Scholes模型。然后，准确地知道基础密度，并由p（s，τ）=σBS给出√τφ日志（秒）- （rd- 射频-σBS）τσBS√τs、 对于s>0，τ>0，（2.17），其中φ（x）是标准正态分布随机变量的密度函数。（2.17）中的密度函数p（s，τ）满足PDEτp=sσBSsp-s（（rd- rf）sp），对于s，τ>0，p（s，0）=δ（s-S） 。该PDE的形式为（2.8），具有空间域的自然限制。请注意，对于数值实验，我们没有应用第1节中的对数变换来获得非常数系数，这使得问题更具挑战性。首先，空间域被截断为[Smin，Smax]=[0，30S]，我们通过使用统一网格的sinh变换定义非统一网格Smin=s<·····<Smax，参见[14]。图1显示了S=100、m=50和S=0到S=5的值的空间网格，以说明点S=S周围较小的网格宽度。应用第2.2小节中的FV离散化，然后得出精确值P（si，τ）的近似值Pi（τ）。s0 100 200 300 400 500图1：Black–Scholes示例中s=Sf周围的非均匀网格图，实际值s=100，m=50。当试图确定数值方法相对于参考解的性能时，重要的是要注意计算值的计算环境，并了解对比较的影响。我们的计算采用64位IEEE浮点算法。由于前向Kolmogorov方程的解是一个概率密度函数，因此解的大小在计算域上变化很大。

这对于初始条件（Dirac delta）尤其如此，对于具有可达到边界的自然有界随机过程，通常情况下，对于0<α的c.f.theSLV模型（1.1），情况也是如此≤ 1/2。由于IEEE floating point有固定长度尾数，因此密度函数无法在整个域内统一表示为高绝对精度。在密度函数较大的地区，只能达到较高的相对精度（正确的位数）。此外，由于数值解是通过隐式时间步进获得的，因此在精确解较小的区域，我们不应期望数值密度具有较高的相对精度。这是因为在求解线性系统时，我们将许多不同量级的项组合在一起，然后将它们相加，这在IEEE算法中会导致相对精度的损失。因此，在比较数值解和参考解时，我们采用了一种混合绝对相对误差度量：当参考解大于1时，我们使用相对误差；如果参考解小于1，则使用绝对误差；我们取整个域的最大误差值。更准确地说，让i（米）=p（si，T）-Pi（T）p（si，T）如果p（si，T）>1，| p（si，T）- Pi（T）|其他。总混合空间误差由以下公式确定：（m） =最大值1≤我≤m级i（m）。1的值有点任意。然而，只要交叉值不太小，结果对交叉值就不太敏感。在实际实验中，通过采用具有大量步长的Crank–Nicolson时间步进方案来近似Pi（T）值，使得时间离散误差可以忽略不计。在图2的左图中，显示了相关情况下的总混合空间误差，其中rd=0.03，rf=0.01，σBS=0.2，T=1，空间网格点的数量m={50，100，…，1000}。

相应的数值解形式=200如图2的右图所示。收敛图清楚地表明，对于当前的初边值问题，theFV离散化是二阶收敛的。作为第二个例子，我们考虑Cox–Ingersoll–Ross（CIR）过程，参见[3]，dVτ=κ（η- Vτ）dτ+ξpVτdWτ，其中κ、η、ξ是严格的正参数。相应的密度函数由，seee给出。g、 [3]，p（v，τ）=ce-u-u（uu）q/2Iq（2√uu），（2.18），其中c=2κξ（1-e-κτ），u=cVe-κτ，u=cv，q=2κηξ- 1,1/m10-310-210-1总体混合空间误差10-710-610-510-4 FV模式的收敛性0 50 100 150 200PN00.0050.010.0150.02数值密度图2：1D Black-Scholes模型内的收敛结果。参数值arerd=0.03，rf=0.01，σBS=0.2，T=1。v0 0.1 0.2图3：CIR示例中v=0和v=v周围的非均匀网格图，实际值v=0.0625，m=50。Iq（·）是第一类q阶的修正贝塞尔函数。请注意，q的值与所谓的费勒条件直接相关，即Vτ=0是可附加的。密度函数（2.18）满足正向Kolmogorov方程τp=vξvp-v（κ（η- v） p），（2.19）对于v，τ>0，其中p（v，0）=δ（v- 五） 。很容易看出，如果伐木条件被违反，即如果q<0，则（2.18）中的密度在v=0时没有定义，并且随着v趋于零，密度函数趋于一致。此外，在v=0左右，PDE（2.19）是强对流占优的，这对数值离散化方法来说非常具有挑战性。域被截断为[Vmin，Vmax]=[0，15]，非均匀网格0=v<v<··<vm=Vmaxis，通过使用均匀底层网格的sinh变换定义。

空间网格如图3所示，V=0.0625、m=50和V=0到V=0.2的值用于净化V=0和V=V周围的较小网格宽度。之后应用FV离散化，从而得出近似值Pj（T），（1≤ j≤ m） 。回想一下，如果q<0，则密度函数趋向于完整，因为v趋向于零，而在v=0时，没有定义精确的密度函数。通过增加空间网格点m的数量，第二个网格点v的值结束为零，并且很难充分比较p（v，T）和p（T）之间的差异。鉴于此，我们选择在相似的空间域上计算误差。如果空间网格点m的总数为50，则设Vlow为最小的非零网格点，即点v。然后，我们用Maxj定义总的混合空间误差≤j≤m级j（m），其中，对于给定的m，jis是最低指数，使得vj≥ 沃兰德j（米）=p（vj，T）-Pj（T）p（vj，T）如果p（vj，T）>1，| p（vj，T）- Pj（T）|其他。1/m10-310-210-1总体混合空间误差10-510-410-310-210-1 FV模式的收敛0.2 0.4 0.6 0.8 1PN012345数值密度图4：CIR模型内的收敛结果。参数由集合A给出。近似值Pj（T）是通过考虑大量时间步长来确定的，因此时间离散化误差可以忽略不计。请注意，确定Vlow的选择m=50并不重要。只要通过实验中考虑的最粗网格之一定义VLOW，数值实验的结论基本不变。对于实际实验，我们考虑了两组参数：κηξVTSet A 5 0.16 0.9 0.0625 0.25 set B 1.15 0.0348 0.39 0.0348 0.25这些参数集取自[7]，也用于[26]。对于集合A，q=0.98，方差过程保持严格正。对于集合B，我们有q=-可达到0.47和Vτ=0。

在图4和图5的左图中，分别显示了集合A和集合B的参数以及空间网格点的数量m={50，100，…，1000}的总混合空间误差。在右图中，m=200时显示了相应的数值解。收敛图表明，对于当前的初边值问题，FV离散化是收敛的。其他实验表明，如果满足Feller条件，FV离散化是二阶收敛的。如果q<0，收敛阶可以降到一阶。此外，所有实验都证实，即使严重违反伐木条件，总数值质量（2.6）始终保持等于1。2.4二维前向Kolmogorov方程在本小节中，一维情况下的FV离散化用于确定一般二维前向Kolmogorov方程的空间离散化（2.1）。假设空间域被截断为[xmin，xmax]×[ymin，ymax]，xd方向上的空间网格点（分别为y方向）由xmin=x<x<…<xm=xmax，分别为ymin=y<y<…<ym=ymax。允许xi=xi- xi-1和yj=yj- yj公司-1be空间网格宽度，其中x=xm+1=y=ym+1=0，确定体积Ohmi、 j：=[xi-0.5，xi+0.5]×[yj-0.5，yj+0.5]，其中xi±0.5，yj±0.5的定义与一维情况类似。1/m10-310-210-1总体混合空间误差10-510-410-310-210-1 FV模式的收敛0.2 0.4 0.6 0.8 1PN0102030数值密度图5：CIR模型内的收敛结果。参数由集合B给出。我们现在可以考虑由pi，j（τ）定义的单元平均值、体积平均值、pi，j（τ）的二维等效值=|Ohmi、 j | ZOhmi、 jp（x，y，τ）dxdy，带|Ohmi、 j |=（xi+0.5）-xi-0.5）（yj+0.5-yj公司-0.5）相应体积的面积。

volumeaverage pi，j（τ）再次是p（xi，yj，τ）的二阶近似值，前提是底层网格是平滑的，并且是通过FV离散化近似的量。已经证实|Ohmi、 j|τpi，j（τ）=Zyj+0.5yj-0.5x个σp- upx=xi+0.5x=xi-0.5dy（2.20a）+Zxi+0.5xi-0.5yσp- upy=yj+0.5y=yj-0.5dx（2.20b）+h[ρσσp]| x=xi+0.5x=xi-0.5iy=yj+0.5y=yj-0.5，（2.20c），二维离散化基于近似值|Ohmi、 j|τpi，j（τ）≈x个σp- upx=xi+0.5x=xi-0.5y=yjyj+yj+1（2.21a）+yσp- upy=yj+0.5y=yj-0.5x=xixi+xi+1（2.21b）+h[ρσσp]| x=xi+0.5x=xi-0.5iy=yj+0.5y=yj-0.5。（2.21c）方程式（2.20）包括几个与第2.1和2.2小节中的flux项类似的flux项。它反映了一个事实，即p在体积上的总积分仅因体积边界上的流体差异而变化。这与（2.3）的一维解释完全类似。如果空间域边界上的总flux为零，即if0=Zymaxyminx个σp- upx=xmaxx=xmindy+Zxmaxxminyσp- upy=ymaxy=ymindx+[ρσσp]| x=xmaxx=xminy=ymaxy=ymin，那么整个区域上p的总积分在时间上是常数。设Pi，j（τ）表示精确值Pi，j（τ）的近似值，用P（τ）表示m×mmatrix，条目为Pi，j（τ）和letP（τ）=vec[P（τ）]，其中vec[·]表示通过将任何给定矩阵的后续列置于彼此下方而将其转化为向量的运算符。粗体符号的引入只是为了表示近似的矩阵形式和向量形式之间的细微差别。

与第2.1和2.2小节中的一维离散化类似，我们通过引入数值fluxespi，j（τ）=[fi-0.5，j（P，τ）- fi+0.5，j（P，τ）]xi+xi+1（2.22a）+[fi，j-0.5（P，τ）- fi，j+0.5（P，τ）]yj+yj+1（2.22b）+Xi，j=0(-1） i+jfm，i-0.5+i，j-0.5+jxi+xi+1yj+yj+1，（2.22c）用于1≤ 我≤ m、 1个≤ j≤ m、 为了便于表述，表示u1，i±0.5，j=u（xi±0.5，yj，τ），σ1，i，j=σ（xi，yj，τ），u2，i，j±0.5=u（xi，yj±0.5，τ），σ2，i，j=σ（xi，yj，τ），σ1，i±0.5，j±0.5=σ（xi±0.5，yj±0.5，τ），σ2，i±0.5，j±0.5=σ（xi±0.5，yj±0.5，τ）。由于（2.21）中流体的实际形式与第2.4节中流体的形式完全相似，我们将数值流体定义为fi±0.5，j（P，τ）=fa，i±0.5，j（P，τ）+fd，i±0.5，j（P，τ），带fa，i-0.5，j（P，τ）=u1，i-0.5，jPi-1，j（τ）+π，j（τ）≈ u1，i-0.5，jp（xi-0.5，yj，τ），和fd，i-0.5，j（P，τ）=σ1，i-1，jPi-1，j（τ）-σ1，i，jPi，j（τ）xi≈ -x个σ（x，yj，τ）p（x，yj，τ）|x=xi-0.5，对于2≤ 我≤ m、 1个≤ j≤ m。