pdeSolver

pdeSolver(coeff, initial, domain, tGrid, [boundary], [method], [setting])

使用有限差分法（FDM）求解一维或二维线性抛物型偏微分方程（PDE）。支持均匀或非均匀结构化网格、Dirichlet/Neumann/Robin 边界条件、一维 obstacle 约束、一维 theta 格式、二维显式格式和 ADI 格式。

可用于以下场景：

• 金融衍生品定价：基于 Black-Scholes 或局部波动率（Local Vol）模型的期权定价，以及美式期权定价。

• 物联网与工程仿真：一维热传导、扩散和对流扩散等。

一维情况下，求解如下方程：

二维情况下，求解如下方程：

求解上述方程时，可配合以下辅助函数生成空间网格、构造求解区域和方程系数，并查看相关配置信息：

构造 domain 和 coeff 后，将其与初始条件 initial、时间网格 tGrid 一并传入 pdeSolver，并按需通过 boundary 指定边界条件，即可进行求解。通过 method 选择数值方法和计算策略，通过 setting 控制输出形式、插值和诊断。

注意： 当前版本不支持 3D PDE、FEM、FVM、非结构网格、强非线性 PDE、二维 obstacle、完整二维隐式格式、定义域外插值和任意非网格时间点的时间插值。

coeff 字典，描述 PDE 方程及其系数。一维情况下由 pdeCoeff1D 构造，二维情况下由 pdeCoeff2D 构造。系数字段的值可以是有限数值标量或用户自定义函数。函数型系数必须在每个求值点返回有限数值标量，且 mass 必须始终为正数。

initial 数值标量、数值向量、数值矩阵或函数，表示 tGrid[0] 时刻的初始条件。在金融场景中通常表示 payoff：

• 一维情况下，数值标量会广播到全部空间节点；数值向量的长度必须等于 x 网格节点数；函数签名为 def(x)。

• 二维情况下，数值标量会广播到全部空间节点；数值矩阵的逻辑形状必须为 nx × ny；函数签名为 def(x, y)。不支持用一维向量描述二维初值。

若使用 Dirichlet 边界，initial 在相应边界节点上的值必须与 tGrid[0] 时刻的边界值一致。二维 Neumann 或 Robin 边界上的节点由求解器根据边界关系回填。

domain 字典，描述空间网格和求解区域。一维情况下由 pdeDomain1D 构造，二维情况下由 pdeDomain2D 构造。

tGrid 一维数值向量，表示时间推进网格。至少包含两个时间点；所有元素必须为有限的非 NULL 值，并严格递增。

boundary 可选参数，数值标量或包含 2 个（一维）或 4 个（二维）元素的向量，表示边界条件。向量中的每个元素可以是数值标量、函数或字典。默认值为 NULL，即所有边界均采用值为 0 的 Dirichlet 边界。传入数值标量时，为所有边界指定相同的常量 Dirichlet 值。

一维情况下，可传入按 [left, right] 排列的二元向量。每个元素可以是：

• 数值标量，表示常量 Dirichlet 边界；

• 签名为 def(t) 的函数，表示随时间变化的 Dirichlet 边界；

• 显式描述 Dirichlet、Neumann 或 Robin 边界的字典。

二维情况下，可传入按 [xMin, xMax, yMin, yMax] 排列的四元向量。xMin/xMax 边界函数的签名为 def(y, t)，yMin/yMax 边界函数的签名为 def(x, t)。

边界字典的字段如下：

边界字典中的值可以是数值标量或符合相应边界函数签名的函数。二维两条 Dirichlet 边在角点相交时，边界值必须在相对容差 1e-12 内一致；一条 Dirichlet 边与 Neumann 或 Robin 边相交时，角点采用 Dirichlet 值。

method 可选参数，STRING 或 SYMBOL 键的字典，表示数值方法和计算策略。默认值为 NULL。不允许包含未定义字段，也不能包含 setting 的字段。

通用字段如下：

一维情况下还支持以下字段：

一维显式格式不能与 method 中的 obstacle 约束或 Rannacher smoothing 组合；PSOR 字段和 penalty 字段不能交叉使用。

二维情况下还支持以下字段：

采用二维显式格式时，method 中不能设置 adiScheme、adiTheta、crossTermTreatment、linearSolver、numThreads 字段，也不能将 rannacherSteps 设为非零值。二维不支持 obstacle 约束；parallel=false 时不能显式设置 numThreads。

setting 可选参数，STRING 或 SYMBOL 键的字典，只控制输出形式、插值和诊断。默认值为 NULL。不允许包含未定义字段，也不能包含 method 的字段。

返回值由 setting 中的 resultMode、output、interpAt 和 calcGreeks 字段决定。

当 resultMode="value"（默认值），且没有指定 interpAt、calcGreeks=false 时，直接返回 PDE 数值解：

二维矩阵的行对应 x 网格，列对应 y 网格，采用 x-fastest:index=i+nx*j 布局。

在 resultMode="value" 下请求插值或 Greeks 时，返回一个字典，按需包含以下字段：

当 resultMode="detail" 时返回字典。始终包含 value；指定 interpAt 时增加 interpValue 和 interpAt；指定 calcGreeks=true 时增加 greeks 和 greekNames；output 为 "all" 或 "selected" 时增加 outputTimes。

当同时指定 diagnostics=true 时，字典还包含：

当 output="final" 且只有一个插值点时，interpValue 为 DOUBLE 标量；多个插值点时为 DOUBLE 向量；当 output="all" 或 output="selected" 时，interpValue 按时间和插值点组成向量或矩阵。

以下期权定价示例以剩余期限 $\tau$ 为时间变量，$\tau=0$ 对应到期时刻，initial 为到期收益。tGrid 从 0 推进至 finalTime，得到估值时点的期权价值。时间以年为单位，利率、分红率和波动率均采用年化值，利率采用连续复利。

输出表中的 pdePrice 为 PDE 数值解；解析解、二叉树或蒙特卡洛结果用于对照。以下输出值作适当舍入。

例1. Black-Scholes 欧式看涨期权。

对不支付分红的标的资产计算欧式看涨期权价格。当前价格和执行价均为 100，无风险利率为 3%，波动率为 20%，剩余期限为 1 年。到期收益为 $\max(S-K,0)$，仅在到期时行权。

空间网格覆盖 $S\in[0,400]$，在执行价附近加密。左端边界为 0，右端边界为 $400-100e^{-0.03\tau}$。通过 setting 的 interpAt 返回当前价格 100 处的插值结果，并与 Black-Scholes 解析解比较。

spot = 100.0
strike = 100.0
rate = 0.03
volatility = 0.20
finalTime = 1
spotMax = 400.0
spotGrid = pdeGrid1D(xMin=0.0, xMax=spotMax, n=801, gridType="sinh", focus=strike, density=3.0)
tGrid = (0..5000) * finalTime / 5000.0
coeff = pdeCoeff1D(
    diff = {S, tau -> 0.02 * S * S },
    conv = {S, tau -> 0.03 * S },
    react = - rate
)
payoff = {S -> max(S - 100.0, 0.0) }
boundary = [
    {tau -> 0.0 },
    {tau -> 400.0 - 100.0 * exp(-0.03 * tau) }
]
setting = dict(keyType=STRING, valueType=ANY)
setting["interpAt"] = spot
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=spotGrid),
    tGrid=tGrid,
    boundary=boundary,
    setting=setting
)
d1 = (log(spot / strike) + (rate + 0.5 * volatility * volatility) * finalTime) / (volatility * sqrt(finalTime))
d2 = d1 - volatility * sqrt(finalTime)
reference = spot * cdfNormal(mean=0.0, stdev=1.0, X=d1) - strike * exp(- rate * finalTime) * cdfNormal(mean=0.0, stdev=1.0, X=d2)
priceError = abs(result[`interpValue] - reference)/reference

table(
    [result["interpValue"]] as pdePrice,
    [reference] as referencePrice,
    [priceError] as relativeError
)

referencePrice 为解析解，relativeError 为数值解相对于解析解的绝对相对误差。

例2. Black-Scholes 美式看跌期权。

计算可在到期前提前行权的美式看跌期权价格。当前价格和执行价均为 100，无风险利率为 5%，波动率为 20%，剩余期限为 1 年，不考虑分红。提前行权收益为 $\max(K-S,0)$，期权价值在各时间层均不能低于该收益。

空间范围为 $[0,300]$，左右边界分别为 100 和 0。将收益函数设为 method 的 obstacle，使用默认的 PSOR 算法处理提前行权约束，并通过 rannacherSteps 平滑初始收益的折点。以 2000 步 CRR 二叉树结果作为数值对照。

spot = 100.0
strike = 100.0
rate = 0.05
volatility = 0.20
finalTime = 1.0
spotGrid = pdeGrid1D(xMin=0.0, xMax=300.0, n=401, gridType="sinh", focus=strike, density=3.0)
tGrid = (0..200) * finalTime / 200.0
coeff = pdeCoeff1D(
    diff={S, tau -> 0.02 * S * S },
    conv={S, tau -> 0.05 * S },
    react=-rate
)
payoff = {S -> max(100.0 - S, 0.0) }
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method["rannacherSteps"] = 2
method["obstacle"] = payoff
method["psorTolerance"] = 1e-10
method["psorMaxIterations"] = 20000
setting["interpAt"] = spot
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=spotGrid),
    tGrid=tGrid,
    boundary=[100.0, 0.0],
    method=method,
    setting=setting
)
def CRR_tree(spot,strike,rate,volatility,finalTime,N){
    dt = finalTime/N
    u = exp(volatility * sqrt(dt))
    d = exp(- volatility * sqrt(dt))
    p = (exp(rate * dt)- d)/(u-d)
    if(p <= 0.0 || p >= 1)
        throw("Invalid risk-neutral probability.")
    discount = exp(-rate*dt)
    j = 0..N
    S = spot * pow(u,j) * pow(d,N-j)
    payoff = strike - S
    V = iif(cond=payoff > 0.0, trueResult=payoff, falseResult=0.0)
    for(i in(N - 1)..0){
        j = 0..i
        S = spot * pow(u,j) * pow(d,i-j)
        continuation = discount * (p * V[j+1] + (1.0 - p) * V[j])
        payoff = strike - S
        exercise = iif(cond=payoff >0.0, trueResult=payoff, falseResult=0.0)
        V = iif(cond=exercise > continuation, trueResult=exercise, falseResult=continuation)
    }
    return V[0]
}
crrReference = CRR_tree(spot=spot, strike=strike, rate=rate, volatility=volatility, finalTime=finalTime, N=2000)
priceError = abs((result["interpValue"] - crrReference)/crrReference)

table(
    [result["interpValue"]] as pdePrice,
    [crrReference] as crrPrice,
    [priceError] as relativeDifference
)

crrPrice 为 CRR 二叉树价格，relativeDifference 为两种数值方法的绝对相对差值。

例3. 局部波动率模型下的欧式看涨期权。

计算执行价为 100、剩余期限为 1 年的欧式看涨期权在标的价格 80、100 和 120 处的价值。无风险利率为 3%，连续分红率为 1%，局部波动率为 $\sigma_{\mathrm{loc}}(S,\tau)=0.2[1+0.2\tanh((S-100)/100)](1+0.1\tau)$。

扩散项随价格和时间变化，对流项为 $(r-q)S=0.02S$。空间范围为 $[0,400]$，左端边界为 0，右端边界为 $400e^{-0.01\tau}-100e^{-0.03\tau}$。分别使用 201、301 和 501 个空间节点及 101、151 和 251 个时间点，比较网格加密后的价格；同时比较同一中等网格上的两种系数预计算策略。

def localVolatility(S, tau) {
    return 0.2 * (1.0 + 0.2 * tanh((S - 100.0) / 100.0)) * (1.0 + 0.1 * tau)
}
def localVolDiffusion(S, tau) {
    sigma = localVolatility(S=S, tau=tau)
    return 0.5 * sigma * sigma * S * S
}
spots = 80.0 100.0 120.0
spotMax = 400.0
coeff = pdeCoeff1D(
    diff=localVolDiffusion,
    conv=def(S, tau) { return 0.02 * S },
    react=-0.03
)
payoff = def(S) { return max(S - 100.0, 0.0) }
boundary = [def(tau) { return 0.0 }, def(tau) { return 400.0 * exp(-0.01 * tau) - 100.0 * exp(-0.03 * tau) }]
coarseGrid = pdeGrid1D(xMin=0.0, xMax=spotMax, n=201, gridType="sinh", focus=100.0, density=3.0)
mediumGrid = pdeGrid1D(xMin=0.0, xMax=spotMax, n=301, gridType="sinh", focus=100.0, density=3.0)
referenceGrid = pdeGrid1D(xMin=0.0, xMax=spotMax, n=501, gridType="sinh", focus=100.0, density=3.0)
coarseMethod = dict(keyType=STRING, valueType=ANY)
coarseMethod["precompute"] = "all"
allMethod = dict(keyType=STRING, valueType=ANY)
allMethod["precompute"] = "all"
sliceMethod = dict(keyType=STRING, valueType=ANY)
sliceMethod["precompute"] = "timeSlice"
referenceMethod = dict(keyType=STRING, valueType=ANY)
referenceMethod["precompute"] = "timeSlice"
setting = dict(keyType=STRING, valueType=ANY)
setting["interpAt"] = spots
coarse = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=coarseGrid),
    tGrid=(0..100) / 100.0,
    boundary=boundary,
    method=coarseMethod,
    setting=setting
)
mediumAll = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=mediumGrid),
    tGrid=(0..150) / 150.0,
    boundary=boundary,
    method=allMethod,
    setting=setting
)
mediumSlice = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=mediumGrid),
    tGrid=(0..150) / 150.0,
    boundary=boundary,
    method=sliceMethod,
    setting=setting
)
reference = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=referenceGrid),
    tGrid=(0..250) / 250.0,
    boundary=boundary,
    method=referenceMethod,
    setting=setting
)
coarseReferenceDifference = max(abs(coarse["interpValue"] - reference["interpValue"]))
mediumReferenceDifference = max(abs(mediumAll["interpValue"] - reference["interpValue"]))
precomputeDifference = max(abs(mediumAll["interpValue"] - mediumSlice["interpValue"]))
prices = mediumAll["interpValue"]

table(
    spots as spot,
    coarse["interpValue"] as coarsePrice,
    prices as mediumPrice,
    reference["interpValue"] as finePrice,
    abs(prices - reference["interpValue"]) as mediumDifference,
    abs(prices - mediumSlice["interpValue"]) as cacheDifference
)

coarsePrice、mediumPrice 和 finePrice 分别对应粗、中、细网格。mediumDifference 为中等网格与细网格价格的绝对差值；cacheDifference 为中等网格上 "all" 与 "timeSlice" 两种预计算策略所得价格的绝对差值。

本例中，中等网格相对细网格的最大差值约为 0.00046813，小于粗网格的约 0.00138244；两种预计算策略的结果一致。

例4. 无分红的向下敲出欧式看涨期权。

计算无返还、连续监测障碍的向下敲出看涨期权。当前标的价格和执行价均为 100，下障碍为 80，无风险利率为 3%，波动率为 20%，剩余期限为 2 年，不考虑分红。存续期间一旦价格触及或跌破 80，期权立即失效且不支付返还；未触及障碍时，到期收益为 $\max(S-K,0)$。

将求解区域限制在 $[80,400]$，在障碍处施加值为 0 的 Dirichlet 边界，右端采用 $400-100e^{-0.03\tau}$。计算当前价格 100 处的期权价值，并与无返还向下敲出看涨期权的解析解比较。

spot = 100.0
strike = 100.0
rate = 0.03
volatility = 0.20
finalTime = 2
spotMax = 400.0
B = 80
spotGrid = pdeGrid1D(xMin=B, xMax=spotMax, n=401, gridType="sinh", focus=strike, density=3.0)
tGrid = (0..200) * finalTime / 200.0
coeff = pdeCoeff1D(
    diff = {S, tau -> 0.02 * S * S },
    conv = {S, tau -> 0.03 * S },
    react = - rate
)
payoff = {S -> max(S - 100.0,0) }
boundary = [
    {tau -> 0.0 },
    {tau -> 400.0 - 100.0 * exp(-0.03 * tau) }
]
setting = dict(keyType=STRING, valueType=ANY)
setting["interpAt"] = spot
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=spotGrid),
    tGrid=tGrid,
    boundary=boundary,
    setting=setting
)
d1 = (log(spot / strike) + (rate + 0.5 * pow(volatility,2)) * finalTime) / (volatility * sqrt(finalTime))
d2 = d1 - volatility * sqrt(finalTime)
d3 = (log((B * B / spot) / strike) + (rate + 0.5 * pow(volatility,2)) * finalTime) / (volatility * sqrt(finalTime))
d4 = d3 - volatility * sqrt(finalTime)
C1 = spot * cdfNormal(mean=0.0, stdev=1.0, X=d1) - strike *exp(- rate * finalTime) * cdfNormal(mean=0.0, stdev=1.0, X=d2)
C2 = B * B / spot * cdfNormal(mean=0.0, stdev=1.0, X=d3) - strike *exp(- rate * finalTime) * cdfNormal(mean=0.0, stdev=1.0, X=d4)
reference = C1 - pow(spot / B ,1 - 2 * rate / (volatility * volatility)) * C2
priceError = abs(result[`interpValue] - reference) / reference

table(
    [result["interpValue"]] as pdePrice,
    [reference] as referencePrice,
    [priceError] as relativeError
)

例5. 现金或无二元看涨期权（Digital）。

计算到期按条件支付固定现金的二元看涨期权。当前标的价格和执行价均为 100，无风险利率为 3%，波动率为 20%，剩余期限为 1 年，不考虑分红。若到期价格严格大于 100，则支付 5；否则支付 0。

空间范围为 $[0,400]$，左端边界为 0，右端边界为 $5e^{-0.03\tau}$。使用 Crank-Nicolson 格式，并设置 rannacherSteps 为 1，以缓解不连续收益引起的数值振荡。与现金或无二元看涨期权的解析解比较。

spot = 100.0
strike = 100.0
rate = 0.03
volatility = 0.20
finalTime = 1
spotMax = 400.0
spotGrid = pdeGrid1D(xMin=0.0, xMax=spotMax, n=401, gridType="sinh", focus=strike, density=3.0)
tGrid = (0..100) * finalTime / 100.0
cashPayoff = 5
coeff = pdeCoeff1D(
    diff = {S, tau -> 0.02 * S * S },
    conv = {S, tau -> 0.03 * S },
    react = - rate
)
payoff = {S -> iif(cond=S > 100.0, trueResult=5, falseResult=0) }
boundary = [
    0.0,
    {tau -> 5.0 * exp(-0.03 * tau) }
]
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method['rannacherSteps'] = 1
method['theta'] = 0.5
setting["interpAt"] = spot
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=spotGrid),
    tGrid=tGrid,
    boundary=boundary,
    method=method,
    setting=setting
)
d2 = (log(spot / strike) + (rate - 0.5 * pow(volatility,2)) * finalTime) / (volatility * sqrt(finalTime))
reference = cashPayoff * exp(- rate * finalTime) * cdfNormal(mean=0.0, stdev=1.0, X=d2)
priceError = abs(result[`interpValue] - reference)/reference

table(
    [result["interpValue"]] as pdePrice,
    [reference] as referencePrice,
    [priceError] as relativeError
)

本例在给定网格下的相对误差约为 0.5742%。收益在执行价处不连续，需通过空间网格加密检查离散误差。

例6. 含连续分红的单边障碍看涨期权（Single Barrier）。

分别计算连续监测、无返还的向下敲出和向上敲出看涨期权。两种情况的当前标的价格和执行价均为 100，无风险利率为 3%，连续分红率为 1%，波动率为 20%，剩余期限为 1 年；对流项系数均为 $(r-q)S$。

（1）向下敲出：下障碍为 80，价格触及或跌破障碍即失效。求解区域为 $[80,400]$，左端边界为 0，右端边界为 $400e^{-0.01\tau}-100e^{-0.03\tau}$。

spot = 100.0
strike = 100.0
barrier = 80.0
rate = 0.03
dividend = 0.01
volatility = 0.20
finalTime = 1.0
spotMax = 400.0
spotGrid = pdeGrid1D(xMin=barrier, xMax=spotMax, n=401, gridType="sinh", focus=strike, density=3.0)
tGrid = (0..100) * finalTime / 100.0
coeff = pdeCoeff1D(
    diff = {S, tau -> 0.02 * S * S },
    conv = {S, tau -> 0.02 * S },
    react = - rate
)
payoff = {S -> max(S - 100.0, 0.0) }
boundary = [
    {tau -> 0.0 },
    {tau -> 400.0 * exp(-0.01 * tau) - 100.0 * exp(-0.03 * tau)}
]
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method["rannacherSteps"] = 2
setting["interpAt"] = spot
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=spotGrid),
    tGrid=tGrid,
    boundary=boundary,
    method=method,
    setting=setting
)
rootTime = volatility * sqrt(finalTime)
mu = (rate - dividend - 0.5 * volatility * volatility) / (volatility * volatility)
x1 = log(spot / strike) / rootTime + (1.0 + mu) * rootTime
y1 = log(barrier * barrier / (spot * strike)) / rootTime + (1.0 + mu) * rootTime
reference = (
    spot * exp(-dividend * finalTime) * (
        cdfNormal(mean=0.0, stdev=1.0, X=x1) -
        pow(barrier / spot, 2.0 * (mu + 1.0)) * cdfNormal(mean=0.0, stdev=1.0, X=y1)) -
    strike * exp(-rate * finalTime) * (
        cdfNormal(mean=0.0, stdev=1.0, X=x1 - rootTime) -
        pow(barrier / spot, 2.0 * mu) * cdfNormal(mean=0.0, stdev=1.0, X=y1 - rootTime))
)
priceError = (abs(result["interpValue"] - reference)) / reference

table(
    [result["interpValue"]] as pdePrice,
    [reference] as referencePrice,
    [priceError] as relativeError
)

（2）向上敲出：上障碍为 130，价格触及或超过障碍即失效。求解区域为 $[0,130]$，两端均采用零 Dirichlet 边界。初值在上障碍处也设为 0，以与边界条件保持一致。两个子例均与相应的解析解比较。

spot = 100.0
strike = 100.0
barrier = 130.0
rate = 0.03
dividend = 0.01
volatility = 0.20
finalTime = 1.0
spotGrid = pdeGrid1D(xMin=0, xMax=barrier, n=401, gridType='sinh', focus=strike, density=3.0)
tGrid = (0..100) * finalTime / 100
coeff = pdeCoeff1D(
    diff = {S,tau -> 0.02 * S * S},
    conv = {S,tau -> 0.02 * S},
    react = -0.03
)
payoff = {S -> iif(cond=S >= 130.0, trueResult=0.0, falseResult=max(S - 100.0, 0))}
boundary = [
    0.0,
    0.0
]
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method['rannacherSteps'] = 2
setting['interpAt'] = spot
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=spotGrid),
    tGrid=tGrid,
    boundary=boundary,
    method=method,
    setting=setting
)
rootTime = volatility * sqrt(finalTime)
mu = (rate - dividend - 0.5 * volatility * volatility) / (volatility * volatility)
x1 = log(spot / strike) / rootTime + (1.0 + mu) * rootTime
x2 = log(spot / barrier) / rootTime + (1.0 + mu) * rootTime
y1 = log(barrier * barrier / (spot * strike)) / rootTime + (1.0 + mu) * rootTime
y2 = log(barrier / spot) / rootTime + (1.0 + mu) * rootTime
reference = (
    exp(-dividend * finalTime) * spot * (
        cdfNormal(mean=0.0, stdev=1.0, X=x1) -
        cdfNormal(mean=0.0, stdev=1.0, X=x2) -
        pow(barrier / spot, 2 * (mu + 1.0))*(cdfNormal(mean=0.0, stdev=1.0, X=y1)-cdfNormal(mean=0.0, stdev=1.0, X=y2))) -
    exp(-rate * finalTime) * strike * (
        cdfNormal(mean=0.0, stdev=1.0, X=x1 - rootTime) -
        cdfNormal(mean=0.0, stdev=1.0, X=x2 - rootTime) -
        pow(barrier / spot, 2 * mu) * (
        cdfNormal(mean=0.0, stdev=1.0, X=y1 - rootTime)-cdfNormal(mean=0.0, stdev=1.0, X=y2 - rootTime)))
)
priceError = abs(result['interpValue'] - reference) / reference

table(
    [result["interpValue"]] as pdePrice,
    [reference] as referencePrice,
    [priceError] as relativeError
)

例7. 触碰后到期支付的 One Touch 期权。

计算向上触碰、到期支付型 One Touch 期权。当前标的价格为 100，上障碍为 125，无风险利率为 3%，波动率为 20%，剩余期限为 1 年，不考虑分红。存续期间若价格触及 125，则在到期时支付 5；始终未触及则支付 0。

求解区域为 $[0,125]$。障碍以下的到期收益为 0，在障碍点的到期收益为 5；左端边界为 0，障碍处边界为 $5e^{-0.03\tau}$，表示触碰后已确定的到期现金支付在当前时刻的现值。与解析解比较。

spot = 100.0
rate = 0.03
volatility = 0.20
finalTime = 1
barrier = 125.0
cashPayoff = 5.0
spotGrid = pdeGrid1D(xMin=0.0, xMax=barrier, n=126, gridType="sinh", focus=spot, density=3.0)
tGrid = (0..100) * finalTime / 100.0
coeff = pdeCoeff1D(
    diff = {S, tau -> 0.02 * S * S },
    conv = {S, tau -> 0.03 * S },
    react = - rate
)
payoff = {S -> iif(cond=S < 125.0, trueResult=0.0, falseResult=5.0)}
boundary = [
    0.0,
    {tau -> 5 * exp(-0.03 * tau) }
]
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method["rannacherSteps"] = 2
setting["interpAt"] = spot
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=spotGrid),
    tGrid=tGrid,
    boundary=boundary,
    method=method,
    setting=setting
)
a  = (rate - 0.5 * pow(volatility,2)) * finalTime
b = log(barrier / spot)
c = volatility * sqrt(finalTime)
d1 = (a - b) / c
d2= (a + b) / c
reference = (
    cashPayoff * exp(- rate * finalTime)*(cdfNormal(mean=0.0, stdev=1.0, X=d1) +
        pow(barrier / spot, 2 * rate / pow(volatility,2)-1) * cdfNormal(mean=0.0, stdev=1.0, X=- d2))
)
priceError = abs(result[`interpValue] - reference) / reference

table(
    [result["interpValue"]] as pdePrice,
    [reference] as referencePrice,
    [priceError] as relativeError
)

例8. 连续算术平均亚式看涨期权。

计算从估值时点开始平均的固定执行价亚式看涨期权。当前标的价格和执行价均为 100，无风险利率为 5%，波动率为 20%，平均期和剩余期限均为 1 年，不考虑分红。到期收益为 $\max(A_T/T-K,0)$，其中 $A_t=\int_0^t S_s\,ds$ 为累计价格积分，本例初始值 $A_0=0$。

通过 $x=(A-KT)/S$ 和 $u(S,A,\tau)=(S/T)v(x,\tau)$ 将问题化为一维 PDE：$v_\tau=\tfrac12\sigma^2x^2v_{xx}+(1-rx)v_x$。在 $[-5,5]$ 上求解，初值为 $\max(x,0)$，左端边界为 0，右端 Neumann 边界为 $v_x=e^{-r\tau}$。在 $x_0=-1$ 处插值后，将结果乘以 $S_0/T$ 还原期权价格。

使用 200000 条蒙特卡洛路径、每条路径 2000 个时间步，以梯形法累计价格积分，作为独立数值对照。固定随机种子以便复现，并输出蒙特卡洛标准误。

spot = 100.0
strike = 100.0
rate = 0.05
volatility = 0.20
finalTime = 1.0
A0 = 0.0
x0 = (A0 - strike * finalTime) / spot
xMin = -5.0
xMax = 5.0
xGrid = pdeGrid1D(xMin=xMin, xMax=xMax, n=201, gridType="sinh", focus=0.0, density=3.0)
tGrid = (0..100) * finalTime / 100.0
coeff = pdeCoeff1D(
    diff = {x, tau -> 0.02 * x * x},
    conv = {x, tau -> 1.0 - 0.05 * x},
    react = 0.0
)
payoff = {x -> max(x, 0.0)}
boundary = [
    0.0,
    dict(keyObj=["type", "value"], valueObj=["neumann", {tau -> exp(- 0.05 * tau)}])
]
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method["rannacherSteps"] = 2
setting["interpAt"] = x0
method["precompute"] = "timeSlice"
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain1D(xGrid=xGrid),
    tGrid=tGrid,
    boundary=boundary,
    method=method,
    setting=setting
)
v0 = result["interpValue"]
price = spot / finalTime * v0
setRandomSeed(seed=20260917)
pathNum = 200000
timeStep = 2000
dt = finalTime / timeStep
integralS = array(dataType=DOUBLE, initialSize=pathNum, capacity=pathNum, defaultValue=0.0)
S = array(dataType=DOUBLE, initialSize=pathNum, capacity=pathNum, defaultValue=spot)
for(i in 1..timeStep){
    z = normal(mean=0.0, std=1.0, count=pathNum)
    SNew = S * exp((rate - 0.5 * volatility * volatility) * dt + volatility * sqrt(dt) * z)
    integralS += 0.5 * (S + SNew) * dt
    S = SNew
}
averageS = integralS / finalTime
payoffMC = max(averageS - strike,0.0)
discount = exp(- rate * finalTime)
priceMC = discount * avg(payoffMC)
mcStdError = discount * std(payoffMC) / sqrt(pathNum)
priceError = abs(price - priceMC) / priceMC

table(
    [price] as pdePrice,
    [priceMC] as mcPrice,
    [mcStdError] as mcStdError,
    [abs(price-priceMC)] as absoluteDifference
)

mcPrice 为蒙特卡洛估计值，mcStdError 为其采样标准误，absoluteDifference 为其与 PDE 价格的绝对差值。本例差值约为 1.46 倍标准误。

例9. Heston 随机波动率模型下的欧式看涨期权。

计算标的价格和瞬时方差共同驱动的欧式看涨期权。当前价格和执行价均为 100，无风险利率为 5%，剩余期限为 1 年，不考虑分红。初始方差 $v_0=0.04$，长期方差 $\theta=0.04$，均值回复速度 $\kappa=2$，方差波动率 $\xi=0.3$，价格与方差的相关系数 $\rho=-0.7$。

代码以 $x=\ln S$ 和方差 $v$ 为二维坐标，价格范围为 $[1,2000]$，方差范围为 $[0,1]$。到期收益为 $\max(e^x-100,0)$；价格方向的左右边界分别为 0 和 $2000-100e^{-0.05\tau}$，方差方向两端采用零 Neumann 近似边界。

使用 MCS ADI 格式处理包含交叉导数的二维 PDE，并在 [log(100), 0.04] 处插值。与相同参数下通过 Heston 特征函数数值积分得到的参考价格 10.394218565150163 比较。

logGrid = (0..200) * ((log(2000.0) - log(1.0)) / 200.0) + log(1.0)
varianceUnit = (0..80) / 80.0
varianceGrid = varianceUnit * varianceUnit
coeff = pdeCoeff2D(
    xx={x, v, t -> 0.5 * v },
    yy={x, v, t -> 0.5 * 0.3 * 0.3 * v },
    xy={x, v, t -> -0.7 * 0.3 * v },
    x={x, v, t -> 0.05 - 0.5 * v },
    y={x, v, t -> 2.0 * (0.04 - v) },
    u=-0.05
)
boundary = [
    0.0,
    def(v, tau) { return 2000.0 - 100.0 * exp(-0.05 * tau) },
    dict(keyObj=["type", "value"], valueObj=["neumann", 0.0]),
    dict(keyObj=["type", "value"], valueObj=["neumann", 0.0])
]
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method["adiScheme"] = "MCS"
method["parallel"] = false
setting["interpAt"] = [log(100.0), 0.04]
method["rannacherSteps"] = 2
method["precompute"] = "timeSlice"
result = pdeSolver(
    coeff=coeff,
    initial=def(x, v) { return max(exp(x) - 100.0, 0.0) },
    domain=pdeDomain2D(xGrid=logGrid, yGrid=varianceGrid),
    tGrid=(0..160) / 160.0,
    boundary=boundary,
    method=method,
    setting=setting
)
reference = 10.394218565150163
price = result["interpValue"]
priceError = (abs(price - reference) / reference)

table(
    [price] as pdePrice,
    [reference] as referencePrice,
    [priceError] as relativeError
)

例10. 双资产篮子看涨期权。

计算两项相关资产等权组成的欧式篮子看涨期权。两项资产当前价格均为 100，权重均为 0.5，执行价为 100，无风险利率为 5%，剩余期限为 1 年，不考虑分红。两项资产的波动率分别为 20% 和 25%，相关系数为 0.3，到期收益为 $\max(0.5S_1+0.5S_2-100,0)$。

两项资产的价格范围均为 $[0,400]$。当某项资产价格为 0 时，边界值由以另一项资产加权价格为标的的 Black-Scholes 看涨期权价格给出；两条价格上界采用外法向导数为 0.5 的 Neumann 条件。二维系数中的 xy 描述资产相关性产生的交叉导数项，使用 MCS ADI 格式求解。

在两项资产价格均为 100 处插值，并与 1000000 组相关终值模拟得到的蒙特卡洛价格比较。固定随机种子并报告采样标准误。

spot1 = 100.0
spot2 = 100.0
strike = 100.0
w1 = 0.5
w2 = 0.5
rate = 0.05
sigma1 = 0.20
sigma2 = 0.25
rho = 0.30
finalTime = 1.0
S1Max = 400.0
S2Max = 400.0
S1Grid = pdeGrid1D(xMin=0.0, xMax=S1Max, n=401, gridType="sinh", focus=spot1, density=3.0)
S2Grid = pdeGrid1D(xMin=0.0, xMax=S2Max, n=401, gridType="sinh", focus=spot2, density=3.0)
tGrid = (0..200) * finalTime / 200.0
coeff = pdeCoeff2D(
    xx = {S1, S2, tau -> 0.5 * 0.20 * 0.20 * S1 * S1},
    yy = {S1, S2, tau -> 0.5 * 0.25 * 0.25 * S2 * S2},
    xy = {S1, S2, tau -> 0.3 * 0.20 * 0.25 * S1 * S2},
    x = {S1, S2, tau -> 0.05 * S1},
    y = {S1, S2, tau -> 0.05 * S2},
    u = - 0.05
)
payoff = {S1, S2 -> max(0.50 * S1 + 0.50 * S2 - 100.0, 0.0)}
def bsCall(S, K, r, sigma, tau) {
    if(tau <= 0.0) {
        return max(S - K, 0.0)
    }
    if(S <= 0.0) {
        return 0.0
    }
    d1 = (log(S / K) + (r + 0.5 * sigma * sigma) * tau) / (sigma * sqrt(tau))
    d2 = d1 - sigma * sqrt(tau)
    return S * cdfNormal(mean=0.0, stdev=1.0, X=d1) - K * exp(-r * tau) * cdfNormal(mean=0.0, stdev=1.0, X=d2)
}
boundary = [
    def(S2, tau) { return bsCall(S=0.5 * S2, K=100.0, r=0.05, sigma=0.25, tau=tau)},
    dict(keyObj=["type", "value"], valueObj=["neumann", 0.5]),
    def(S1, tau) {return bsCall(S=0.5 * S1, K=100.0, r=0.05, sigma=0.20, tau=tau)},
    dict(keyObj=["type", "value"], valueObj=["neumann", 0.5])
]
method = dict(keyType=STRING, valueType=ANY)
setting = dict(keyType=STRING, valueType=ANY)
method["adiScheme"] = "MCS"
method["parallel"] = false
setting["interpAt"] = [spot1, spot2]
method["rannacherSteps"] = 2
method["precompute"] = "timeSlice"
result = pdeSolver(
    coeff=coeff,
    initial=payoff,
    domain=pdeDomain2D(xGrid=S1Grid, yGrid=S2Grid),
    tGrid=tGrid,
    boundary=boundary,
    method=method,
    setting=setting
)
price = result["interpValue"]
setRandomSeed(seed=20260917)
N = 1000000
z1 = normal(mean=0.0, std=1.0, count=N)
zIndependent = normal(mean=0.0, std=1.0, count=N)
z2 = rho * z1 + sqrt(1.0 - rho * rho) * zIndependent
S1T = spot1 * exp((rate - 0.5 * sigma1 * sigma1) * finalTime + sigma1 * sqrt(finalTime) * z1)
S2T = spot2 * exp((rate - 0.5 * sigma2 * sigma2) * finalTime + sigma2 * sqrt(finalTime) * z2)
payoffMC = max(w1 * S1T + w2 * S2T - strike,0.0)
discount = exp(-rate * finalTime)
priceMC = discount * avg(payoffMC)
relativeError = abs(result["interpValue"] - priceMC) / priceMC

mcStdError = discount * std(payoffMC) / sqrt(N)

table(
    [price] as pdePrice,
    [priceMC] as mcPrice,
    [mcStdError] as mcStdError,
    [abs(price-priceMC)] as absoluteDifference
)

absoluteDifference 为 PDE 价格与蒙特卡洛估计值的绝对差值，本例约为 0.31 倍蒙特卡洛标准误。

例11. 求解一维热方程 $u_t=0.2u_{xx}$。空间范围为 $[0,1]$，左右两端采用零 Dirichlet 边界，初值为 $u(x,0)=x(1-x)$。

heatXGrid = (0..40) / 40.0
heatTGrid = (0..800) * (0.05 / 800.0)
heatCoeff = pdeCoeff1D(diff=0.2)
heatInitial = heatXGrid * (1.0 - heatXGrid)
heatDomain = pdeDomain1D(xGrid=heatXGrid)

heatResult = pdeSolver(
    coeff=heatCoeff,
    initial=heatInitial,
    domain=heatDomain,
    tGrid=heatTGrid,
    boundary=0.0
)

heatResult

例12. 求解二维热方程 $u_t=0.05(u_{xx}+u_{yy})$。区域为 $[0,1]\times[0,1]$，四条边均采用零 Dirichlet 边界。

plateXGrid = (0..10) / 10.0
plateYGrid = (0..10) / 10.0
plateTGrid = (0..100) / 1000.0
plateCoeff = pdeCoeff2D(xx=0.05, yy=0.05)

def plateInitial(x, y) {
    return x * (1.0 - x) * y * (1.0 - y)
}

plateDomain = pdeDomain2D(xGrid=plateXGrid, yGrid=plateYGrid)
plateResult = pdeSolver(
    coeff=plateCoeff,
    initial=plateInitial,
    domain=plateDomain,
    tGrid=plateTGrid,
    boundary=0.0
)

plateResult

例13. 求解二维泊松型问题。边界向量按 xMin、xMax、yMin、yMax 排列，其中 xMin 和 yMin 使用 Robin 边界，xMax 和 yMax 使用 Neumann 边界。

boundaryXGrid = 0.0 0.2 0.8 1.5
boundaryYGrid = 0.0 0.3 1.0 1.8
boundaryTGrid = 0.0 0.05 0.1 0.15 0.2

boundaryDomain = pdeDomain2D(xGrid=boundaryXGrid, yGrid=boundaryYGrid)
boundaryCoeff = pdeCoeff2D(xx=1.0, yy=1.0, source=-6.0)
def boundaryInitial(x, y) {
    return x * x + 2.0 * y * y + 3.0 * x + 4.0 * y + 5.0
}

xMinRobin = dict(STRING, ANY)
xMinRobin["type"] = "robin"
xMinRobin["alpha"] = 2.0
xMinRobin["beta"] = 0.5
xMinRobin["value"] = def(y, t) {
    return 2.0 * (2.0 * y * y + 4.0 * y + 5.0) - 1.5
}

xMaxNeumann = dict(STRING, ANY)
xMaxNeumann["type"] = "neumann"
xMaxNeumann["value"] = 6.0

yMinRobin = dict(STRING, ANY)
yMinRobin["type"] = "robin"
yMinRobin["alpha"] = 1.5
yMinRobin["beta"] = 0.75
yMinRobin["value"] = def(x, t) {
    return 1.5 * (x * x + 3.0 * x + 5.0) - 3.0
}

yMaxNeumann = dict(STRING, ANY)
yMaxNeumann["type"] = "neumann"
yMaxNeumann["value"] = 11.2

boundaryConditions = [xMinRobin, xMaxNeumann, yMinRobin, yMaxNeumann]

boundaryResult = pdeSolver(
    coeff=boundaryCoeff,
    initial=boundaryInitial,
    domain=boundaryDomain,
    tGrid=boundaryTGrid,
    boundary=boundaryConditions
)

boundaryResult

例14. 使用支持 JIT 的 DolphinDB Server，在 method 中将 jit 字段设为 true，为系数、初值和边界回调请求 JIT 编译。先定义具名用户自定义函数，再将函数名传入相应参数。

jitXGrid = (0..40) / 40.0
jitTGrid = (0..200) * (0.05 / 200.0)

def jitDiff(x, t) {
    return 0.2
}

def jitInitial(x) {
    return x * (1.0 - x)
}

def jitBoundary(t) {
    return 0.0
}

jitCoeff = pdeCoeff1D(diff=jitDiff)
jitDomain = pdeDomain1D(xGrid=jitXGrid)
jitMethod = dict(keyType=STRING, valueType=ANY)
jitMethod["jit"] = true

jitResult = pdeSolver(
    coeff=jitCoeff,
    initial=jitInitial,
    domain=jitDomain,
    tGrid=jitTGrid,
    boundary=[jitBoundary, jitBoundary],
    method=jitMethod
)

jitResult

例15. 设置随时间变化的 Dirichlet 边界。求解 $\frac{\partial u}{\partial t}=0.2\frac{\partial^2u}{\partial x^2}+1$，空间范围为 [0,1]，初值为 $u(x,0)=x$，左右边界分别为 $u(0,t)=t$ 和 $u(1,t)=1+t$。左端直接传入 def(t) 函数，右端通过字典设置，展示两种等价的描述方式。初值与 t=0 时的边界值一致，解析解为 $u(x,t)=x+t$。

rampXGrid = (0..4) / 4.0
rampTGrid = 0.0 0.05 0.1
rampDomain = pdeDomain1D(xGrid=rampXGrid)
rampCoeff = pdeCoeff1D(diff=0.2, source=1.0)

def rampLeftValue(t) {
    return t
}
def rampRightValue(t) {
    return 1.0 + t
}
rampRight = dict(keyType=STRING, valueType=ANY)
rampRight["type"] = "dirichlet"
rampRight["value"] = rampRightValue

rampResult = pdeSolver(
    coeff=rampCoeff,
    initial=rampXGrid,
    domain=rampDomain,
    tGrid=rampTGrid,
    boundary=[rampLeftValue, rampRight]
)
round(rampResult, 6)
// output: [0.1,0.35,0.6,0.85,1.1]

例16. 设置非零 Neumann 边界，并组合 Dirichlet 与 Neumann 边界。求解 $\frac{\partial u}{\partial t}=0.2\frac{\partial^2u}{\partial x^2}-0.4$，空间范围为 [0,1.5]，初值为 $u(x,0)=x^2+3x+5$，解析解不随时间变化。Neumann 的 value 指定外法向导数：左端为 $-\left.\frac{\partial u}{\partial x}\right|_{x=0}=-3$，右端为 $\left.\frac{\partial u}{\partial x}\right|_{x=1.5}=6$。以下使用非均匀网格；默认的二阶边界差分要求至少 4 个节点。第二次求解将左端改为固定值 5，右端仍使用 Neumann 边界。

gradientXGrid = 0.0 0.2 0.8 1.5
gradientTGrid = 0.0 0.05 0.1
gradientInitial = gradientXGrid * gradientXGrid + 3.0 * gradientXGrid + 5.0
gradientDomain = pdeDomain1D(xGrid=gradientXGrid)
gradientCoeff = pdeCoeff1D(diff=0.2, source=-0.4)

gradientLeft = dict(keyType=STRING, valueType=ANY)
gradientLeft["type"] = "neumann"
gradientLeft["value"] = -3.0
gradientRight = dict(keyType=STRING, valueType=ANY)
gradientRight["type"] = "neumann"
gradientRight["value"] = 6.0

gradientResult = pdeSolver(
    coeff=gradientCoeff,
    initial=gradientInitial,
    domain=gradientDomain,
    tGrid=gradientTGrid,
    boundary=[gradientLeft, gradientRight]
)
round(gradientResult, 6)
// output: [5,5.64,8.04,11.75]

mixedResult = pdeSolver(
    coeff=gradientCoeff,
    initial=gradientInitial,
    domain=gradientDomain,
    tGrid=gradientTGrid,
    boundary=[5.0, gradientRight]
)
round(mixedResult, 6)
// output: [5,5.64,8.04,11.75]

例17. 用函数设置 Robin 边界字典中的 alpha、beta 和 value。求解 $\frac{\partial u}{\partial t}=0.2\frac{\partial^2u}{\partial x^2}+1$，初值为 $u(x,0)=x$，解析解为 $u(x,t)=x+t$。两端均取 $\alpha(t)=1+t$、$\beta(t)=0.5+0.1t$。按照 $\alpha(t)u+\beta(t)\frac{\partial u}{\partial n}=\mathrm{value}(t)$，左端的 $\mathrm{value}(t)=(1+t)t-(0.5+0.1t)$，右端的 $\mathrm{value}(t)=(1+t)^2+(0.5+0.1t)$。减号和加号分别来自左、右端的外法向方向。

robinXGrid = (0..4) / 4.0
robinTGrid = 0.0 0.05 0.1
robinDomain = pdeDomain1D(xGrid=robinXGrid)
robinCoeff = pdeCoeff1D(diff=0.2, source=1.0)

def robinAlpha(t) {
    return 1.0 + t
}
def robinBeta(t) {
    return 0.5 + 0.1 * t
}
def robinLeftValue(t) {
    return (1.0 + t) * t - (0.5 + 0.1 * t)
}
def robinRightValue(t) {
    return (1.0 + t) * (1.0 + t) + (0.5 + 0.1 * t)
}

robinLeft = dict(keyType=STRING, valueType=ANY)
robinLeft["type"] = "robin"
robinLeft["alpha"] = robinAlpha
robinLeft["beta"] = robinBeta
robinLeft["value"] = robinLeftValue
robinRight = dict(keyType=STRING, valueType=ANY)
robinRight["type"] = "robin"
robinRight["alpha"] = robinAlpha
robinRight["beta"] = robinBeta
robinRight["value"] = robinRightValue

robinResult = pdeSolver(
    coeff=robinCoeff,
    initial=robinXGrid,
    domain=robinDomain,
    tGrid=robinTGrid,
    boundary=[robinLeft, robinRight]
)
round(robinResult, 6)
// output: [0.1,0.35,0.6,0.85,1.1]

相关函数： pdeCoeff1D、pdeCoeff2D、pdeDomain1D、pdeDomain2D、pdeGrid1D、pdeInfo