0 引言

用随机微分方程来描述客观现象在实际应用中有着越来越重要的作用^[1-6], 通常用数值方法得到的解来近似方程的解析解.人们主要研究数值求解方法的收敛性^[7]和稳定性, 研究的较多的是Euler方法、Milstein方法及Runge-Kutta方法.近年来又提出Heun方法, 本文在该数值方法的基础上构造出θ-Heun方法, 并研究了该方法用于求解标量自治随机微分方程的收敛性.

1 随机微分方程及数值求解方法 1.1 随机微分方程及性质

一维标量自治的随机微分方程^[8]:

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{d}}X\left( t \right) = f\left( {X\left( t \right)} \right){\rm{d}}t + g\left( {X\left( t \right)} \right){\rm{d}}w\left( t \right), }&{t \in [{t_0}, T], }\\ {X({t_0}) = {X_0}, }&{X \in {\bf{R}}, } \end{array}} \right. $

(1)

其中:f(X), g(X)在[t₀, T]上连续可测; E|X₀|² < ∞; w(t)是标准的winner过程, 且与X₀相互独立.当t>0, 步长h>0时, 增量Δw(t)=w(t+h)-w(t)是一列服从正态分布N(0, h), 且相互独立的随机变量, 有如下性质^[9].

1) E[∫_a^bg(X(t))dw(t)]=0, a < b.

2) E[|∫_a^bg(X(t))dw(t)|²]=∫_c^dE[|g(X(t))|²]dt, c < d.

1.2 θ-Heun方法

定义1  求解随机微分方程(1)的θ-Heun方法为

$ \begin{array}{l} {X_{n + 1}} = {X_n} + h[\left( {1 - \theta } \right)f({X_n}) + \theta f({X_n} + hf({X_n}))] + \\ g({X_n})\Delta {w_n}, \theta \in \left[ {0, 1} \right]. \end{array} $

(2)

2 θ-Heun方法的收敛性

定义2^[8]  记X(t_n)为随机微分方程(1)在网格点t_n处的精确值, X_n为通过方法(2)得到的随机微分方程(1)的数值解, $ \tilde X({t_{n + 1}})$是通过方法(2)在X(t_n)处迭代1步得到方程(1)的数值解, 即$\tilde X({t_{n + 1}}) $=X(t_n)+ϕ(h, X(t_n), Δw_n), 其中ϕ(h, X(t_n), Δw_n)=h[(1-θ)f(X(t_n))+θf(X(t_n)+hf(X(t_n)))]+g(X(t_n))Δw_n, 则θ-Heun方法的局部截断误差和整体误差分别为:

$ \begin{array}{l} {\delta _{n + 1}} = X({t_{n + 1}}) - \tilde X({t_{n + 1}}), \;\;{\varepsilon _{n + 1}} = \\ X({t_{n + 1}}) - {X_{n + 1}}, \;n = 0, 1, \cdots , N - 1. \end{array} $

定义3^[6]  如果存在常数C>0(C与h无关), 则$\mathop {\max }\limits_{0 \le n \le N - 1} |E({\delta _{n + 1}})| \le C{h^{{p_1}}} $, $\mathop {\max }\limits_{0 \le n \le N - 1} {[E|{\delta _{n + 1}}{|^2}]^{\frac{1}{2}}} \le C{h^{{p_2}}} $, $\mathop {\max }\limits_{0 \le n \le N - 1} {[E|{\varepsilon _{n + 1}}{|^2}]^{\frac{1}{2}}} \le C{h^p} $, p₁、p₂分别是数值方法在均值意义、均方意义上的收敛阶, p是均方强收敛阶.

引理1(Hölder不等式)^[6]  如果函数f(x)与g(x)在区间[a, b]上连续可积, 则

$ \smallint _a^bf\left( x \right)g\left( x \right){\rm{d}}x \le {(\smallint _a^b{\left( {f\left( x \right)} \right)^p}{\rm{d}}x)^{\frac{1}{p}}}{(\smallint _a^b{\left( {g\left( x \right)} \right)^q}{\rm{d}}x)^{\frac{1}{q}}}. $

引理2^[6]  如果f(X)、g(X)满足以下条件:

1) 线性增长条件.存在一个正的常数L₁, 使

|f(X)|∨|g(X)|≤L₁(1+|X|), 或|f(X)|²∨|g(X)|²≤L₁(1+|X|²).

2) Lipschitz条件.存在一个正的常数L₂, 使|f(X)-f(Y)|∨|g(X)-g(Y)|≤L₂|X-Y|, ∀X, Y∈R. “∨”表示两函数中较大者, 若f(X), g(X)满足条件1)和2), 则方程(1)的解满足

ⅰ)$ E{\left| {X\left( t \right)} \right|^2} \le E\mathop {\sup }\limits_{{t_0} \le t \le T} {\left| {X\left( t \right)} \right|^2} \le {C_1}$.

ⅱ) ∀t₀≤s≤t≤T, E|X(t)-X(s)|²≤C₂(t-s).

其中:C₁、C₂为仅依赖于初值X₀和T的常数.

引理3^[10]  对于方程(1), 当h→0且ε₀=0时, 若$\mathop {\max }\limits_{0 \le n \le N - 1} (E|{\varepsilon _{n + 1}}|) = {\rm O}({h^{{p_3}}}) $, $\mathop {\max }\limits_{0 \le n \le N - 1} (E|{\varepsilon _n}{|^2}) = {\rm O}({h^{{p_4}}}) $, 那么就有p₄=2p₃.

定理1  若f(X)、g(X)满足引理1和2, 对于θ-Heun方法, p₁=2, p₂=1, p=1(步长h≤1).

下面证明p₁=2.

证明

$ \begin{array}{l} {\delta _{n + 1}} = X({t_{n + 1}}) - \tilde X({t_{n + 1}}) = \\ [X({t_n}) + \smallint _{{t_n}}^{{t_{n + 1}}}f\left( {X\left( s \right)} \right){\rm{d}}s + \smallint _{{t_n}}^{{t_{n + 1}}}g\left( {X\left( s \right)} \right){\rm{d}}w] - \\ [X({t_n}) + h[\left( {1 - \theta } \right)f(X({t_n})) + \theta f(X({t_n}) + hf(X({t_n})))] + \\ g(X({t_n}))\Delta {w_n}] = \smallint _{{t_n}}^{{t_{n + 1}}}[f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s + \\ \smallint _{{t_n}}^{{t_{n + 1}}}[g\left( {X\left( s \right)} \right) - g(X({t_n}))]{\rm{d}}w - h\theta (f(X({t_n}) + hf(X({t_n}))) - f(X({t_n}))). \end{array} $

(3)

对上式两端同时取期望有

$ \begin{array}{l} E({\delta _{n + 1}}) \le E\smallint _{{t_n}}^{{t_{n + 1}}}[f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s + E\smallint _{{t_n}}^{{t_{n + 1}}}\\ [g\left( {X\left( s \right)} \right) - g(X({t_n}))]{\rm{d}}w + h\theta E\\ (f(X({t_n}) + hf(X({t_n}))) - f(X({t_n}))). \end{array} $

根据w(t)的性质1)可知

$ \begin{array}{l} E\smallint _{{t_n}}^{{t_{n + 1}}}[g\left( {X\left( s \right)} \right) - g(X({t_n}))]{\rm{d}}w\left( s \right) = 0, \\ E({\delta _{n + 1}}) \le E\smallint _{{t_n}}^{{t_{n + 1}}}[f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s + \\ hE(f(X({t_n}) + hf(X({t_n}))) - f(X({t_n}))), \\ |E({\delta _{n + 1}})| \le |E\smallint _{{t_n}}^{{t_{n + 1}}}[f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s| + \\ h|E(f(X({t_n}) + hf(X({t_n}))) - f(X({t_n})))|. \end{array} $

A)根据引理2的(2)和(3)式及w(t)的性质1)有

$ \begin{array}{l} |E\smallint _{{t_n}}^{{t_{n + 1}}}[f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s| \le {L_2}E\smallint _{{t_n}}^{{t_{n + 1}}}\\ |X\left( s \right) - X({t_n})|{\rm{d}}s = {L_2}E\smallint _{{t_n}}^{{t_{n + 1}}}|X({t_n}) + \smallint _{{t_n}}^s\\ f\left( {X\left( u \right)} \right){\rm{d}}u + \smallint _{{t_n}}^sg\left( {X\left( u \right)} \right){\rm{d}}w\left( u \right) - X({t_n})\\ |{\rm{d}}s \le {L_2}E\smallint _{{t_n}}^{{t_{n + 1}}}\smallint _{{t_n}}^s|f\left( {X\left( u \right)} \right)|{\rm{d}}u{\rm{d}}s + {L_2}E\smallint _{{t_n}}^{{t_{n + 1}}}\smallint _{{t_n}}^s\\ |g\left( {X\left( u \right)} \right)|{\rm{d}}w\left( u \right){\rm{d}}s \le {L_1}{L_2}E\smallint _{{t_n}}^{{t_{n + 1}}}\smallint _{{t_n}}^s\left( {1 + |X\left( u \right)|} \right)\\ {\rm{d}}u{\rm{d}}s \le {L_1}{L_2}E\smallint _{{t_n}}^{{t_{n + 1}}}\smallint _{{t_n}}^s(1 + \frac{1}{2}({1^2} + {\left| {X\left( u \right)} \right|^2})){\rm{d}}u{\rm{d}}s \le \\ {L_1}{L_2}\smallint _{{t_n}}^{{t_{n + 1}}}\smallint _{{t_n}}^s(3 + E{\left| {X\left( u \right)} \right|^2}){\rm{d}}u{\rm{d}}s \le {L_1}{L_2}\smallint _{{t_n}}^{{t_{n + 1}}}\smallint _{{t_n}}^s(3 + {C_1}){\rm{d}}u{\rm{d}}s \le {L_1}{L_2}(3 + {C_1}){h^2}. \end{array} $

B) 根据引理2有

$ \begin{array}{l} h|E(f(X({t_n}) + hf(X({t_n}))) - f(X({t_n})))| \le {L_2}hE|\\ X({t_n}) + hf(X({t_n})) - X({t_n})| \le {L_1}{L_2}{h^2}E(1 + |X({t_n})|)\\ \le {L_1}{L_2}{h^2}E(1 + \frac{1}{2}({1^2} + |X({t_n}){|^2})) \le {L_1}{L_2}{h^2}\\ (3 + E|X({t_n}){|^2}) \le {L_1}{L_2}(3 + {C_1}){h^2}, \end{array} $

从而|E(δ_n+1)|≤2L₁L₂(3+C₁)h², 记C₃=2L₁L₂(3+C₁), 则$\mathop {\max }\limits_{0 \le n \le N - 1} |E({\delta _{n + 1}})| \le {C_3}{h^2} $, 由定义3知, p₁=2.

下面证明p₂=1.

证明  由(3)式及(a+b+c)²≤3a²+3b²+3c²得

$ \begin{array}{l} |{\delta _{n + 1}}{|^2} \le 3|\smallint _{{t_n}}^{{t_{n + 1}}}[f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s{|^2} + 3|\smallint _{{t_n}}^{{t_{n + 1}}}\\ [g\left( {X\left( s \right)} \right) - g(X({t_n}))]{\rm{d}}w\left( s \right){|^2} + 3{h^2}{\theta ^2}|f(X({t_n}) + \\ hf(X({t_n}))) - f(X({t_n})){|^2}, E|{\delta _{n + 1}}{|^2} \le 3E|\smallint _{{t_n}}^{{t_{n + 1}}}\\ [f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s{|^2} + 3E|\smallint _{{t_n}}^{{t_{n + 1}}}[g\left( {X\left( s \right)} \right) - \\ g(X({t_n}))]{\rm{d}}w\left( s \right){|^2} + 3{h^2}E|f(X({t_n}) + hf(X({t_n}))) - f(X({t_n})){|^2}.\\ {\rm{A}})\;E|\smallint _{{t_n}}^{{t_{n + 1}}}[f\left( {X\left( s \right)} \right) - f(X({t_n}))]{\rm{d}}s{|^2} \le E{(\smallint _{{t_n}}^{{t_{n + 1}}}|f\left( {X\left( s \right)} \right) - f(X({t_n}))|{\rm{d}}s)^2}. \end{array} $

(4)

由Hölder不等式得$ \smallint _{{t_n}}^{{t_{n + 1}}}|f\left( {X\left( s \right)} \right) - f(X({t_n}))|{\rm{d}}s \le $$ {(\smallint _{{t_n}}^{{t_{n + 1}}}|f\left( {X\left( s \right)} \right) - f(X({t_n})){|^2}{\rm{d}}s)^{\frac{1}{2}}}{(\smallint _{{t_n}}^{{t_{n + 1}}}{1^2}{\rm{d}}s)^{\frac{1}{2}}}$, 代入(4)式, 再由引理2有

$ \begin{array}{l} h\smallint _{{t_n}}^{{t_{n + 1}}}E|f\left( {X\left( s \right)} \right) - f(X({t_n})){|^2}{\rm{d}}s \le {L_2}^2h\smallint _{{t_n}}^{{t_{n + 1}}}\\ E|X\left( s \right) - X({t_n}){|^2}{\rm{d}}s \le {L_2}^2{C_2}h\smallint _{{t_n}}^{{t_{n + 1}}}(s - {t_n}){\rm{d}}s \le {L_2}^2{C_2}{h^3}. \end{array} $

B) 根据w(t)的性质2)及引理2有

$ \begin{array}{l} E|\smallint _{{t_n}}^{{t_{n + 1}}}[g\left( {X\left( s \right)} \right) - g(X({t_n}))]{\rm{d}}w\left( s \right){|^2} = \smallint _{{t_n}}^{{t_{n + 1}}}\\ E|g\left( {X\left( s \right)} \right) - g(X({t_n})){|^2}{\rm{d}}s \le {L_2}^2\smallint _{{t_n}}^{{t_{n + 1}}}E\\ |X\left( s \right) - X({t_n}){|^2}{\rm{d}}s \le {L_2}^2{C_2}\smallint _{{t_n}}^{{t_{n + 1}}}(s - {t_n}){\rm{d}}s = {L_2}^2{C_2}{h^2}. \end{array} $

C) 由引理2可得

$ \begin{array}{l} E|f[X({t_n}) + hf(X({t_n}))] - f(X({t_n})){|^2} \le {L_2}^2E\\ |X({t_n}) + hf(X({t_n})) - X({t_n}){|^2} = {L_2}^2{h^2}E|f(X({t_n})){|^2}\\ \le {L_2}^2{L_1}{h^2}E(1 + |X({t_n}){|^2}) = {L_2}^2{L_1}{h^2}(1 + E(|X({t_n}){|^2}))\\ \le {L_2}^2{L_1}{h^2}(1 + {C_1}). \end{array} $

由A)、B)、C)得

$ \begin{array}{l} E|{\delta _{n + 1}}{|^2} \le 3{L_2}^2{C_2}{h^3} + 3{L_2}^2{C_2}{h^2} + 3{L_2}^2\\ {L_1}{h^4}(1 + {C_1}) \le 3{L_2}^2(2{C_2} + {L_1}(1 + {C_1})){h^2}.{\rm{ }} \end{array} $

记C₄=3L₂²(C₂+L₁(1+C₁)), 则E|δ_n+1|²≤C₄h², 从而$\mathop {\max }\limits_{0 \le n \le N - 1} {[E|{\delta _{n + 1}}{|^2}]^{\frac{1}{2}}} \le \sqrt {{C_4}} h $.由定义3知, p₂=1.

下面证明p=1.

证明  根据定义1得ε_n+1=X(t_n+1)-X_n+1=[X(t_n+1)-${\tilde X}$(t_n+1)]+[${\tilde X}$(t_n+1)-X_n+1]=δ_n+1+[${\tilde X}$(t_n+1)-X_n+1]=δ_n+1+[X(t_n)+h[(1-θ)f(X(t_n))+θf(X(t_n)+hf(X(t_n)))]+g(X(t_n))Δw_n]-[X_n+h[(1-θ)f(X_n)+θf(X_n+hf(X_n))]+g(X_n)Δw_n]=δ_n+1+ε_n+(g(X(t_n))-g(X_n))Δw_n+h(1-θ)(f(X(t_n))-f(X_n))+hθ(f(X(t_n)+hf(X(t_n)))-f(X_n+hf(X_n)))E|ε_n+1|≤E|ε_n|+E|δ_n+1|+h(1-θ)E|f(X(t_n))-f(X_n)|+E|(g(X(t_n))-g(X_n))Δw_n|+hθE|f(X(t_n)+hf(X(t_n)))-f(X_n+hf(X_n))|.

A) 由引理2的2)及ε_n的定义得h(1-θ)E|f(X(t_n))-f(X_n)|≤L₂hE|X(t_n)-X_n|=L₂hE|ε_n|.

B) 由引理2的2)及w(t)与X_n的独立性^[8]得

$ E|(g(X({t_n})) - g({X_n}))\Delta {w_n}| \le {L_2}E|(X({t_n}) - {X_n})\Delta {w_n}| = {L_2}E|{\varepsilon _n}\left| E \right|\Delta {w_n}|. $

因为Δw(t)~N(0, h), 根据正态分布的定义得$E|\Delta {w_n}| = \sqrt {\left( {2h/\pi } \right)} < 1 $, 则

$ E|(g(X({t_n})) - g({X_n}))\Delta {w_n}| \le {L_2}E|{\varepsilon _n}|. $

C) 由引理2的2)及ε_n的定义得

$ \begin{array}{l} h\theta E|f(X({t_n}) + hf(X({t_n}))) - f({X_n} + hf({X_n}))| \le h{L_2}E\\ |X({t_n}) + hf(X({t_n})) - ({X_n} + hf({X_n}))| = h{L_2}E|X({t_n}) - \\ {X_n} + h(f(X({t_n})) - f({X_n}))| \le h{L_2}E|{\varepsilon _n}| + {h^2}{L_2}E|\\ f(X({t_n})) - f({X_n})| \le h{L_2}E|{\varepsilon _n}| + {h^2}{L_2}^2E|{\varepsilon _n}| = \\ (h{L_2} + {h^2}{L_2}^2)E|{\varepsilon _n}\left| {, E} \right|{\varepsilon _{n + 1}}| \le (1 + {L_2} + (2{L_2} + {L_2}^2h)h)\\ E|{\varepsilon _n}\left| { + E} \right|{\delta _{n + 1}}| \le (1 + {L_2} + (2{L_2} + {L_2}^2)h)E|{\varepsilon _n}\left| { + E} \right|{\delta _{n + 1}}|. \end{array} $

记C₅=2L₂+L₂², 显然C₅>0且是与h无关的常数, 则E|ε_n+1|≤(1+L₂+C₅h)E|ε_n|+E|δ_n+1|.

反复应用上述关系可得

$ \begin{array}{l} E|{\varepsilon _{n + 1}}| \le {(1 + {L_2} + {C_5}h)^{n + 1}}E|{\varepsilon _0}| + \\ (\mathop {\max }\limits_{0 \le n \le N} E[|{\delta _{n + 1}}|\left] ) \right[{(1 + {L_2} + {C_5}h)^n} - 1]/({L_2} + {C_5}h).{\rm{ }} \end{array} $

由引理3, 当ε₀=0时, E|ε₀|=0, 再由(1)的证明上式变为

$ \begin{array}{l} E|{\varepsilon _{n + 1}}| = (\mathop {\max }\limits_{0 \le n \le N} E|{\delta _{n + 1}}|)[{(1 + {L_2} + {C_5}h)^n} - 1]/\\ ({L_2} + {C_5}h) \le ({C_3}{h^2}/({L_2} + {C_5}))h[{(1 + {L_2} + {C_5}h)^n} - 1]\\ \le h\frac{{{C_3}}}{{{L_2} + {C_5}}}({{\rm{e}}^{{C_5}T}} - 1) = O\left( h \right). \end{array} $

由引理3得p₃=1, 故p₄=2, 即存在不依赖于h的常数C₆, 使得$ \mathop {\max }\limits_{{t_0} \le n \le T} {(E|{\varepsilon _n}{|^2})^{\frac{1}{2}}} \le {C_6}h$, 从而p=1.

3 数值实验

对于实验方程

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{d}}X\left( t \right) = \lambda X{\rm{d}}t + \mu X{\rm{d}}w\left( t \right), }&{t \in \left[ {0, T} \right].}\\ {X\left( 0 \right) = {X_0} = 1, }&{\lambda , \mu \in {\bf{R}}.} \end{array}} \right. $

取$ \lambda = - 5, \mu = \frac{1}{5}, \theta = \frac{1}{5}, {X_0} = 1, T = 1$, 在[0, 1]区间上, 选取3个不同的步长h分别为1、1/2、1/4, 利用θ-Heun方法和Heun方法进行迭代, 运用Matlab结果如图 1和图 2所示.从图 1与图 2的对比可以得出:用θ-Heun方法解上述方程的数值解更逼近精确解, θ-Heun方法的平均误差为6.919 1×10^-6, Heun方法的平均误差为1.228 2×10^-4, 验证了由θ-Heun方法得到的数值解的误差与Heun方法相比较小.