pdeSolver

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

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

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

$M(x,t)\frac{\partial u}{\partial t}=A(x,t)\frac{\partial^2u}{\partial x^2}+B(x,t)\frac{\partial u}{\partial x}+C(x,t)u+D(x,t)$

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

$M(x,y,t)\frac{\partial u}{\partial t}=A(x,y,t)\frac{\partial^2u}{\partial x^2}+B(x,y,t)\frac{\partial^2u}{\partial y^2}+C(x,y,t)\frac{\partial^2u}{\partial x\partial y}+D(x,y,t)\frac{\partial u}{\partial x}+E(x,y,t)\frac{\partial u}{\partial y}+F(x,y,t)u+G(x,y,t)$

注意： 当前版本不支持 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 可选参数，表示边界条件。默认值为 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 的字段。

通用字段如下：

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

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

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

二维显式格式不能设置 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 时，字典还包含：

在 final 输出且只有一个插值点时，interpValue 为 DOUBLE 标量；多个插值点时为 DOUBLE 向量；all/selected 输出时按时间和插值点组成向量或矩阵。

例1. 求解一维热方程 $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

例2. 求解二维热方程 $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

例3. 求解二维泊松型问题。边界向量按 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

例4. 在 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

例5. 设置随时间变化的 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]

例6. 设置非零 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]

例7. 用函数设置 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