求解二维变系数G-方程的高精度交替方向隐格式
High Accuracy Alternating Direction Implicit Scheme for Solving Two-Dimensional Variable Coefficient G-Equations
摘要: 本文提出了求解二维变系数G方程的高精度交替方向隐式(ADI)格式,该方法在物理、金融和计算数学中具有重要应用。我们分析了该格式的稳定性,并进行了数值实验,比较了不同格式的绝对误差,验证了ADI格式的有效性。
Abstract: In this paper, we propose high-accuracy Alternating Direction Implicit (ADI) scheme for solving two-dimensional variable coefficient G-equations, which have significant applications in physics, finance and computational mathematics. The stability of the scheme is analyzed, and numerical experiment is conducted to compare absolute errors across different schemes, confirming the effectiveness of the scheme.
文章引用:张宇, 许乾海, 李洋. 求解二维变系数G-方程的高精度交替方向隐格式[J]. 理论数学, 2025, 15(2): 268-276. https://doi.org/10.12677/pm.2025.152067

1. 引言

2007年,Peng [1]通过G-热方程首次提出了G-正态分布和G-期望的概念。随后,Peng [2]基于这一框架推导了次线性期望下的大数定律和中心极限定理,并证明了该框架下的弱收敛性。2010年Peng [3]给出了G-期望和G-方程的联系,Hu, Ji, Peng, Song [4]和Wang [5]进一步证明了G-期望框架下的非线性Feynman-Kac公式,建立了G-布朗运动驱动的随机微分方程与偏微分方程之间的关系。G-热方程(1)作为一种特殊的HJB [6]方程受到了广泛的研究,

U t G( 2 U x 2 )=f( t,x ) (1)

Duncan [7]和Peng [8]证明了G-热方程粘性解的存在性和唯一性,Hu [9]研究了一类特殊G-热方程的解与常微分方程解的关系,并由此得出了该类G-热方程的解。然而,G-热方程仅是G-微分方程中的一种特殊情况,现实中的模型通常比G-热方程复杂得多。目前,对于一般形式的G-方程解析解的计算仍缺乏有效的方法。因此,数值方法成为求解G-方程的常用工具,并被广泛应用。

宋彬彬[10]求解了一维G-方程的Crank-Nicolson格式,Yang和Zhao [11]针对一维非齐次G-热方程应用显式差分格式、半隐式差分格式、隐式差分格式和Crank-Nicolson格式,在对应的数值格式的构造中考虑了迭代系数的改进,王晓莹[12]对传统的有限差分格式进行改进,计算了所改进数值格式的齐次G-热方程的数值解。曹海峰[13]讨论了一维G-热方程的差分方法并给出了相应的数值模拟,还给出了二维G-热方程的显式差分格式、隐式差分格式、半隐式差分格式以及Crank-Nicolson格式,尽管未对这些差分格式进行相关的理论分析和数值模拟,但这些格式的提出为后续研究提供了有价值的出发点。

在许多领域中,二维变系数G-方程扮演着重要的角色,尤其是在物理学和金融学中。这类方程通常用于描述带有空间变异性的现象,例如在热传导、流体力学等物理问题中,二维变系数方程被广泛应用来刻画扩散过程的时空变化。例如,Batchelor [14]介绍了变系数方程在复杂流动模型中的应用,特别是在处理具有空间和时间变异性的流动问题时。在金融领域,二维变系数G-方程同样具有重要的应用,尤其是在期权定价模型中。Black和Scholes [15]提出的经典期权定价模型通过随机微分方程描述资产价格的动态,虽然原始模型假设了常数波动率,但在实际市场中,波动率通常是变化的,因而引入变系数方程来更好地捕捉这种波动性。Heston [16]则提出了一种具有随机波动率的期权定价模型,进一步拓展了这一理论,表明二维变系数G-方程在复杂金融模型中的不可或缺性。

目前,G-方程数值方法的研究大多集中在常系数上。然而,在实际应用中,例如Hu [17]提到的G-期望框架下的Ornstein-Uhlenback过程模型的定价问题,可能会出现变系数。变系数G方程的数值计算方法是一项重要且可推广的工作。我们发现,大多数学者都专注于研究一维变系数G方程式,而二维变系数G-领域仍需要进一步研究。本文将研究二维变系数G-方程的一般形式:

U t a( t,x,y )G( 2 U x 2 )b( t,x,y )G( 2 U y 2 )=f( t,x,y ) (2)

本文的主要工作和创新点如下:提出了一般形式的二维变系数G-方程的ADI格式,对所提差分格式的稳定性进行了严格的理论证明,通过数值实验得到了ADI格式和显式欧拉格式、C-N格式不同节点处的误差。通过计算结果的比较,得到ADI格式误差较小、精度更高。

2. 预备知识(G-期望空间)

本文中作者将研究在G框架下的二维G-方程,对此,我们将沿用文献[1]-[3]中的关于G-期望空间的相关定义以及应用条件。

定义2.1 [1] Ω 是一给定集合, 是定义在 Ω 上的实值函数所组成的一个线性空间,并且满足以下的条件:

a. 每一个实值的常数c都在 中;

b. 如果 X( ) ,则 | X( ) |

那么我们把 中的函数称为随机变量,而称二元组 ( Ω, ) 为随机变量空间。

定义2.2 [2]定义在随机变量空间 上的满足以下性质的泛函 E:R

(I) 单调性:若 AB ,则有 E[ A ]E[ B ]

(II) 保常数性:对任意 γR 都有 E[ γ ]=γ

称三元组 ( Ω,,E ) 为非线性期望空间。称E为一个非线性期望。

E还满足:

(III) 次可加性: E[ A+B ]E[ A ]+E[ B ] A,B

(IV) 正齐次性: E[ ηA ]=ηE[ A ] η0

( Ω,,E ) 为次线性期望空间,E ( Ω, ) 上的次线性期望。

定义2.3 [3] (G-正态分布)

我们称次线性期望空间 ( Ω,,E ) 中的n维随机向量 A=( a 1 ,, a n ) 服从G-正态分布,如果对任意的 α,β>0 ,有:

αA+β A ¯ = d α 2 + β 2 A,

其中 A ¯ A的任意独立复制。

定义2.4 [3] (G-分布)

我们称次线性期望空间 ( Ω,,E ) 中的n维随机向量 ( X,η ) 服从G-分布,如果其满足

( αX+β X ¯ , α 2 η+ β 2 η ¯ ) = d ( α 2 + β 2 X,( α 2 + β 2 )η ),

其中 ( X ¯ , η ¯ ) ( X,η ) 的独立复制。

3. 二维变系数G-方程ADI格式的稳定性分析和数值模拟

本节我们主要研究二维变系数G-方程的ADI数值格式,并对其相应的稳定性进行了分析,在本节最后通过数值算例验证了差分格式的有效性。

3.1. ADI格式

首先我们给出差分方程的ADI格式。

3.1.1. 差分格式的建立

对于一类二维变系数G-方程:

U t a( t,x,y )G( 2 U x 2 )b( t,x,y )G( 2 U y 2 )=f( t,x,y ), (3)

其中 a( t,x,y ) b( t,x,y ) f( t,x,y ) 在定义域上是连续的,并且有 a( t,x,y )>0 b( t,x,y )>0 G( ξ )= 1 2 ( σ ¯ 2 ξ + σ _ 2 ξ ) ξ + =max{ 0,ξ } ξ =min{ 0,ξ }

方程(3)在空间上通常是无界的,而在数值计算过程当中通常是在有界区域中求解,所以我们通过给出第一边值条件对方程进行边界截断,并对G函数进行形式转换,截断后的方程形式如下:

{ t U 1 2 a( t,x,y ) σ 2 ( xx U ) xx U 1 2 b( t,x,y ) σ 2 ( yy U ) yy U=f( t,x,y ), x U| x=l = ϕ 1 ( t ), x U| x=l = ϕ 2 ( t ), y U| y=l = ϕ 3 ( t ), y U| y=l = ϕ 4 ( t ), U| t=0 =ψ( x,y ), (4)

其中:

( t,x,y )[ 0,T ]×[ l,l ]×[ l,l ], σ 2 ( x )={ σ ¯ 2 ,x0 σ _ 2 ,x<0 .

下面主要根据方程(4)给出方程的差分格式。

在求解区间 ( t,x,y )D=[ 0,T ]×[ l,l ]×[ l,l ] 上,将求解区间进行网格化分解,取正整数MN,在 [ l,l ] 上进行M等分,记 h= 2l/M 为空间步长;在 [ 0,T ] 上进行N等分,记 τ=T/N 为时间步长;记 x i =l+ih( 0iM ) y j =l+jh( 0jM ) t n =nτ( 0nN )

交替方向隐式法(ADI)将第n层到n + 1层计算分为两步:先由第n层到n + 1/2层,对 xx U 用向后差分逼近,对 yy U 用向前差分逼近,然后由第n + 1/2层到n + 1层,对 xx U 用向前差分逼近,对 yy U 用向后差分逼近。在处理变系数时,关键是在离散化过程中正确地选择系数的取值。在ADI方法中,通常会在时间步的中间点(如n + 1/2)处取值,以确保数值解的稳定性和准确性。此时,对于定义在D上的网格函数: U i,j n =u( t n , x i , y j ) 。我们使用如下符号:

a n =a( t n , x i , y j ), b n =b( t n , x i , y j ), f i,j n =f( t n , x i , y j )

则G-方程的ADI格式为:

{ δ τ U i,j n+ 1 2 1 2 a n+ 1 2 σ 2 ( δ x 2 U i,j n+ 1 2 ) δ x 2 U i,j n+ 1 2 1 2 b n σ 2 ( δ y 2 U i,j n ) δ y 2 U i,j n = f i,j n , δ τ U i,j n+1 1 2 a n+ 1 2 σ 2 ( δ x 2 U i,j n+ 1 2 ) δ x 2 U i,j n+ 1 2 1 2 b n+1 σ 2 ( δ y 2 U i,j n+1 ) δ y 2 U i,j n+1 = f i,j n+ 1 2 , δ x U 1,j n+1 = ϕ 1 ( t n+1 ), δ x U N x ,j n+1 = ϕ 2 ( t n+1 ), δ y U i,1 n+1 = ϕ 3 ( t n+1 ), δ y U i, N y n+1 = ϕ 4 ( t n+1 ), U i,j 1 =ψ( x i , y j ). (5)

其中 i,j=1,2,,M;n=1,,N 。并且有:

a n+ 1 2 = a n + a n+1 2 , δ τ U i,j n+ 1 2 = U i,j n+ 1 2 U i,j n τ 2 , δ τ U i,j n+1 = U i,j n+1 U i,j n+ 1 2 τ 2 ;

δ x 2 U i,j n+1 = U i+1,j n+1 2 U i,j n+1 + U i1,j n+1 h 2 , δ y 2 U i,j n+1 = U i,j+1 n+1 2 U i,j n+1 + U i,j1 n+1 h 2 ;

δ x U i,j n+1 = U i+1,j n+1 U i1,j n+1 2h , δ y U i,j n+1 = U i,j+1 n+1 U i,j1 n+1 2h .

这样,我们就建立了对于方程(3) ADI差分格式。下面我们对其稳定性进行分析。

3.1.2. ADI差分格式的稳定性

假设 a( t,x,y ):=a( t ) b( t,x,y ):=b( t ) ,仅依赖于t。由极值定理[17],我们可以得到如下定理:

定理3.1对于差分格式(5),我们以 r n+ 1 2 = σ 2 ( δ x 2 u i,j n+ 1 2 )τ/ 2 h 2 r ˜ n = σ 2 ( δ y 2 u i,j n )τ/ 2 h 2 ,定义最大模 u = max 0i,jM | u i,j | f = max 0i,jM | f i,j | 。若差分格式(5)满足条件:

τ h 2 1 σ ¯ 2 max t[ 0,T ] { a( t ),b( t ) }

则有

u n+1 u 0 + k=0 n τ 2 ( f k + f k+ 1 2 ),

即差分格式(5)是稳定的[17]

证明:方程(5)可改写成如下形式:

U i,j n+ 1 2 U i,j n = a n+ 1 2 r n+ 1 2 ( U i+1,j n+ 1 2 2 U i,j n+ 1 2 + U i1,j n+ 1 2 )+ b n r ˜ n ( U i,j+1 n 2 U i,j n + U i,j1 n )+ τ 2 f i,j n , (6)

U i,j n+1 U i,j n+ 1 2 = a n+ 1 2 r n+ 1 2 ( U i+1,j n+ 1 2 2 U i,j n+ 1 2 + U i1,j n+ 1 2 )+ b n+1 r ˜ n+1 ( U i,j+1 n+1 2 U i,j n+1 + U i,j1 n+1 )+ τ 2 f i,j n+ 1 2 , (7)

整理(6) (7)式可得,

( 1+2 a n+ 1 2 r n+ 1 2 ) U i,j n+ 1 2 =( 12 b n r ˜ n ) U i,j n + a n+ 1 2 r n+ 1 2 ( U i+1,j n+ 1 2 + U i1,j n+ 1 2 )+ b n r ˜ n ( U i,j+1 n + U i,j1 n )+ τ 2 f i,j n , (8)

( 1+2 b n+1 r ˜ n+1 ) U i,j n+1 =( 12 a n+ 1 2 r n+ 1 2 ) U i,j n+ 1 2 + a n+ 1 2 r n+ 1 2 ( U i+1,j n+ 1 2 + U i1,j n+ 1 2 )+ b n+1 r ˜ n+1 ( U i,j+1 n+1 + U i,j1 n+1 )+ τ 2 f i,j n+ 1 2 , (9)

当方程满足定理条件时,则有:

12 b n r ˜ n 0,12 a n+ 1 2 r n+ 1 2 0

式子(8) (9)两边同时取绝对值,由绝对值不等式可得

( 1+2 a n+ 1 2 r n+ 1 2 )| U i,j n+ 1 2 |( 12 b n r ˜ n )| U i,j n |+ a n+ 1 2 r n+ 1 2 ( | U i+1,j n+ 1 2 |+| U i1,j n+ 1 2 | )+ b n r ˜ n ( | U i,j+1 n |+| U i,j1 n | )+ τ 2 | f i,j n |, (10)

( 1+2 b n+1 r ˜ n+1 )| U i,j n+1 |( 12 a n+ 1 2 r n+ 1 2 )| U i,j n+ 1 2 |+ a n+ 1 2 r n+ 1 2 ( | U i+1,j n+ 1 2 |+| U i1,j n+ 1 2 | )+ b n+1 r ˜ n+1 ( | U i,j+1 n+1 |+| U i,j1 n+1 | )+ τ 2 | f i,j n+ 1 2 |, (11)

不等式两边取无穷模[17]可得:

( 1+2 a n+ 1 2 r n+ 1 2 ) U n+ 1 2 ( 12 b n r ˜ n ) U n +2 a n+ 1 2 r n+ 1 2 U n+ 1 2 +2 b n r ˜ n U n + τ 2 f n , (12)

( 1+2 b n+1 r ˜ n+1 ) U n+1 ( 12 a n+ 1 2 r n+ 1 2 ) U n+ 1 2 +2 a n+ 1 2 r n+ 1 2 U n+ 1 2 +2 b n+1 r ˜ n+1 U n+1 + τ 2 f n+ 1 2 , (13)

消去多余项后,整理可得

u n+1 u n + τ 2 f n + τ 2 f n+ 1 2

递推可得:

u n+1 u 0 + k=0 n τ 2 ( f k + f k+ 1 2 )

定理得证。

3.2. 数值模拟

算例3.2.1我们考虑下面的二维变系数G-方程:

{ t u 1 2 ( e x + x 2 +t )G( xx u ) 1 2 ( e y + y 2 +t )G( yy u )=f( t,x,y ) u| t=0 =sin( x+y ). (14)

假设方程的精确解为 sin( x+y+t ) ,为使等式成立,此时有:

f=cos( x+y+t ) 1 2 ( e x + x 2 +t ) σ 2 [ sin( x+y+t ) ][ sin( x+y+t ) ] 1 2 ( e y + y 2 +t ) σ 2 [ sin( x+y+t ) ][ sin( x+y+t ) ]

我们令 x,y[ l,l ] t[ 0,T ] σ[ 0.2,1 ] ,其中 l=2π T=2 。首先取 M=40 N=512 进行数值模拟,下面给出此模型的精确解和数值解的比较,图1图2绘制了 t=2 时方程(14)解析解和ADI格式下数值解。

Figure 1. Analytical solution

1. 解析解

Figure 2. Numerical solution of ADI scheme

2. ADI格式数值解

表1给出了部分计算时得到的部分数值结果,数值解很好地逼近解析解,表2比较了ADI格式与显式欧拉格式、C-N格式的绝对误差,结果表明ADI格式的误差低于显式欧拉格式以及C-N格式,表3给出了不同网格精度下的部分网格点ADI格式的绝对误差以及收敛阶,随着网格空间划分的不断精细,误差也逐渐降低,由此得出ADI格式的有效性。

Table 1. Numerical solutions, analytical solutions, and absolute errors between some grid points in ADI scheme at t = 2

1. t = 2时,部分网格点ADI格式的数值解,解析解以及它们之间的绝对误差

u( x i , y j ,1 )

数值解

解析解

绝对误差

u( x 8 , y 8 ,1 )

1.102× 10 1

1.148× 10 1

4.544× 10 3

u( x 16 , y 16 ,1 )

9.235× 10 1

9.093× 10 1

1.417× 10 2

u( x 24 , y 24 ,1 )

6.871× 10 1

6.768× 10 1

1.030× 10 2

u( x 32 , y 32 ,1 )

4.933× 10 1

4.910× 10 1

2.263× 10 3

u( x 40 , y 40 ,1 )

9.810× 10 1

9.802× 10 1

7.318× 10 4

Table 2. Absolute errors of ADI scheme, explicit Euler scheme, and C-N scheme for some grid points at t = 2

2. t = 2时,部分网格点ADI格式、显式欧拉格式和C-N格式的绝对误差

u( x i , y j ,1 )

显式欧拉格式

C-N格式

ADI格式

u( x 8 , y 8 ,1 )

9.210× 10 3

5.051× 10 3

4.544× 10 3

u( x 16 , y 16 ,1 )

5.372× 10 2

2.372× 10 2

1.417× 10 2

u( x 24 , y 24 ,1 )

3.831× 10 2

1.831× 10 2

1.030× 10 2

u( x 32 , y 32 ,1 )

8.356× 10 3

2.588× 10 3

2.263× 10 3

u( x 40 , y 40 ,1 )

5.526× 10 3

1.259× 10 3

7.318× 10 4

Table 3. Absolute error and convergence order of ADI scheme for partial grid points with different grid accuracies at x = 0, y = 0, t = 2

3. x = 0, y = 0, t = 2时,不同网格精度下的部分网格点ADI格式的绝对误差和收敛阶

空间划分

绝对误差

收敛阶

M=20

4.587× 10 3

M=40

1.151× 10 3

1.995

M=60

4.770× 10 4

2.053

M=80

2.689× 10 4

2.054

M=100

1.812× 10 4

2.028

4. 结论

二维变系数G-方程在金融和物理等领域具有重要应用。然而,由于变系数和多维问题的复杂性,理论分析和数值模拟比传统的二维G-方程复杂得多。本文构造了二维变系数G-方程的ADI格式,并进行了详细的稳定性分析。同时,通过数值算例比较了不同方案下不同网格点的数值解和误差,直观地证明了ADI格式的有效性。

基金项目

本论文由上海理工大学教师发展研究项目(编号:CFTD2023YB38)资助。

NOTES

*通讯作者。

参考文献

[1] Peng, S. (2007) G-Expectation, G-Brownian Motion and Related Stochastic Calculus of Itô Type. In: Benth, F.E., et al., Eds., Stochastic Analysis and Applications: The Abel Symposium 2005, Springer, 541-567.
https://doi.org/10.1007/978-3-540-70847-6_25
[2] Peng, S. (2019) Law of Large Numbers and Central Limit Theorem under Nonlinear Expectations. Probability, Uncertainty and Quantitative Risk, 4, Article No. 4.
https://doi.org/10.1186/s41546-019-0038-2
[3] Peng, S. (2019) Nonlinear Expectations and Stochastic Calculus under Uncertainty. Springer.
[4] Hu, M., Ji, S., Peng, S. and Song, Y. (2014) Comparison Theorem, Feynman-Kac Formula and Girsanov Transformation for BSDEs Driven by G-Brownian Motion. Stochastic Processes and Their Applications, 124, 1170-1195.
https://doi.org/10.1016/j.spa.2013.10.009
[5] Peng, S. and Wang, F. (2015) BSDE, Path-Dependent PDE and Nonlinear Feynman-Kac Formula. Science China Mathematics, 59, 19-36.
https://doi.org/10.1007/s11425-015-5086-1
[6] Yong, J. and Zhou, X. (2012) Stochastic Controls: Hamiltonian Systems and HJB Equations. Springer Science & Business Media.
[7] Pasik-Duncan, B. and Duncan, T.E. (2001) Stochastic Controls: Hamiltonian Systems and HJB Equations. IEEE Transactions on Automatic Control, 46, 1846-1846.
https://doi.org/10.1109/tac.2001.964706
[8] Peng, S. (1992) A Generalized Dynamic Programming Principle and Hamilton-Jacobi-Bellman Equation. Stochastics and Stochastic Reports, 38, 119-134.
https://doi.org/10.1080/17442509208833749
[9] Hu, M. (2012) Explicit Solutions of the-Heat Equation for a Class of Initial Conditions. Nonlinear Analysis: Theory, Methods & Applications, 75, 6588-6595.
https://doi.org/10.1016/j.na.2012.08.002
[10] 宋彬彬. G-方程的数值方法[D]: [硕士学位论文]. 苏州: 苏州大学, 2014.
[11] Yang, J. and Zhao, W. (2016) Numerical Simulations for G-Brownian Motion. Frontiers of Mathematics in China, 11, 1625-1643.
https://doi.org/10.1007/s11464-016-0504-9
[12] 王晓莹. 基于有限差分法的G-热方程数值解法[D]: [硕士学位论文]. 济南: 山东大学, 2021.
[13] 曹海峰. HJB偏微分方程的数值计算[D]: [硕士学位论文]. 济南: 山东大学, 2009.
[14] Batchelor, G.K. (2000) An Introduction to Fluid Dynamics. Cambridge University Press.
https://doi.org/10.1017/cbo9780511800955
[15] Black, F. and Scholes, M. (1973) The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81, 637-654.
https://doi.org/10.1086/260062
[16] Heston, S.L. (1993) A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial Studies, 6, 327-343.
https://doi.org/10.1093/rfs/6.2.327
[17] 李荣华, 刘波. 微分方程数值解法[M]. 第4版. 北京: 高等教育出版社, 2009: 83-139.