pdeSolver

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

Solves one-dimensional or two-dimensional linear parabolic partial differential equations (PDEs) with the finite difference method (FDM). It supports uniform and non-uniform structured grids, Dirichlet/Neumann/Robin boundary conditions, 1D obstacle constraints, 1D theta schemes, and 2D explicit and ADI schemes.

Applications include:

• Financial derivatives pricing: option pricing under the Black-Scholes or local volatility (Local Vol) model, and American option pricing.

• IoT and engineering simulations: one-dimensional heat conduction, diffusion, and convection-diffusion.

For a one-dimensional domain, the function solves:

$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)$

For a two-dimensional domain, the function solves:

$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)$

The following helper functions generate spatial grids, construct the domain and coefficient dictionaries, and inspect their configuration for solving these equations:

After constructing domain and coeff, pass them to pdeSolver with the initial condition initial and time grid tGrid, specifying boundary conditions through boundary as needed. Use method to select numerical methods and computational strategies, and setting to control output formats, interpolation, and diagnostics.

Note: The current version does not support 3D PDEs, FEM, FVM, unstructured grids, strongly nonlinear PDEs, 2D obstacles, a fully implicit 2D scheme, extrapolation outside the domain, or interpolation at arbitrary off-grid times.

coeff is a dictionary specifying the PDE and its coefficients. Construct it with pdeCoeff1D for a one-dimensional domain or pdeCoeff2D for a two-dimensional domain. Each coefficient field can contain a finite numeric scalar or a user-defined function. A function-valued coefficient must return a finite numeric scalar at every evaluation point, and mass must remain positive.

initial is a numeric scalar, numeric vector, numeric matrix, or function specifying the initial condition at tGrid[0]. In financial applications, it typically represents the payoff:

• For a 1D domain, a numeric scalar is broadcast to all spatial nodes; a numeric vector must have the same length as the x grid; a function must have the signature def(x).

• For a 2D domain, a numeric scalar is broadcast to all spatial nodes; a numeric matrix must have the logical shape nx × ny; a function must have the signature def(x, y). A one-dimensional vector cannot be used as a 2D initial condition.

For a Dirichlet boundary, initial at the corresponding boundary nodes must match the boundary value at tGrid[0]. Nodes on 2D Neumann or Robin boundaries are filled by the solver from the boundary relations.

domain is a dictionary specifying the spatial grid and solution region. Construct it with pdeDomain1D for a one-dimensional domain or pdeDomain2D for a two-dimensional domain.

tGrid is a one-dimensional numeric vector specifying the time grid. It must contain at least two finite, non-NULL values in strictly increasing order.

boundary (optional) is a numeric scalar or a vector containing 2 elements for a one-dimensional domain or 4 elements for a two-dimensional domain, specifying the boundary conditions. Each vector element can be a numeric scalar, function, or dictionary. The default value is NULL, which applies a zero Dirichlet condition to every boundary. A numeric scalar applies the same constant Dirichlet value to every boundary.

For a 1D domain, specify a two-element vector in [left, right] order. Each element can be:

• A numeric scalar specifying a constant Dirichlet condition;

• A function with the signature def(t) specifying a time-dependent Dirichlet condition;

• A dictionary explicitly specifying a Dirichlet, Neumann, or Robin condition.

For a 2D domain, specify a four-element vector in [xMin, xMax, yMin, yMax] order. Functions for xMin/xMax must have the signature def(y, t), and functions for yMin/yMax must have the signature def(x, t).

A boundary dictionary uses the following fields:

Values in a boundary dictionary can be numeric scalars or functions with the applicable boundary signature. When two Dirichlet boundaries meet at a 2D corner, their values must agree within a relative tolerance of 1e-12. When a Dirichlet boundary meets a Neumann or Robin boundary, the Dirichlet value is used at the corner.

method (optional) is a dictionary with STRING or SYMBOL keys specifying the numerical method and execution policies. The default value is NULL. Unknown fields are not allowed, and fields from setting cannot be placed in method.

The common fields are:

For a 1D domain, the following fields are also supported:

A 1D explicit scheme cannot be combined with an obstacle constraint specified in method or with Rannacher smoothing. PSOR fields and penalty fields cannot be mixed.

For a 2D domain, the following fields are also supported:

For a 2D explicit scheme, method cannot contain the adiScheme, adiTheta, crossTermTreatment, linearSolver, or numThreads fields, and rannacherSteps cannot be nonzero. A 2D obstacle constraint is not supported. When parallel=false, numThreads cannot be specified explicitly.

setting (optional) is a dictionary with STRING or SYMBOL keys controlling only the return format, interpolation, and diagnostics. The default value is NULL. Unknown fields are not allowed, and fields from method cannot be placed in setting.

The return value is determined by the resultMode, output, interpAt, and calcGreeks fields in setting.

When resultMode="value" (the default), interpAt is not specified, and calcGreeks=false, the function directly returns the PDE solution:

Rows of a 2D matrix correspond to the x grid and columns correspond to the y grid, with layout x-fastest:index=i+nx*j.

When interpolation or Greeks are requested in resultMode="value", the function returns a dictionary containing the following fields as applicable:

When resultMode="detail", the function returns a dictionary. It always contains value; adds interpValue and interpAt when interpAt is specified; adds greeks and greekNames when calcGreeks=true; and adds outputTimes when output is "all" or "selected".

When diagnostics=true is also specified, the dictionary additionally contains:

When output="final" and there is a single interpolation point, interpValue is a DOUBLE scalar. For multiple points, it is a DOUBLE vector. When output="all" or output="selected", interpValue is a vector or matrix indexed by time and interpolation point.

The option pricing examples below use time to maturity $\tau$ as the time variable. At $\tau=0$, initial specifies the payoff at maturity. Advancing tGrid from 0 to finalTime gives the option value at the valuation time. Time is measured in years, and interest rates, dividend yields, and volatilities are annualized. Interest rates are continuously compounded.

In the output tables, pdePrice is the PDE price, compared with an analytical, binomial-tree, or Monte Carlo result. Output values below are rounded.

Example 1. Black-Scholes European call option.

Price a European call on a non-dividend-paying asset. The current asset price and strike are both 100, the risk-free rate is 3%, volatility is 20%, and time to maturity is 1 year. The payoff is $\max(S-K,0)$, with exercise permitted only at maturity.

The spatial grid covers $S\in[0,400]$ and is refined near the strike. The left boundary is 0, and the right boundary is $400-100e^{-0.03\tau}$. Set interpAt in setting to obtain the interpolated price at the current asset price of 100, and compare it with the Black-Scholes analytical price.

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 is the analytical price, and relativeError is the absolute relative error against it.

Example 2. Black-Scholes American put option.

Price an American put that permits early exercise. The current asset price and strike are both 100, the risk-free rate is 5%, volatility is 20%, and time to maturity is 1 year, with no dividends. The option value must remain at least as large as the exercise payoff $\max(K-S,0)$ at every time level.

The spatial domain is $[0,300]$, with boundary values of 100 and 0. Set obstacle in method to the payoff function to apply the early-exercise constraint with the default PSOR solver. Use rannacherSteps to smooth the initial payoff kink. Compare the result with a 2000-step CRR binomial tree.

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 is the CRR binomial-tree price, and relativeDifference is the absolute relative difference between the two numerical methods.

Example 3. European call under a local volatility model.

Price a European call with strike 100 and time to maturity of 1 year at asset prices of 80, 100, and 120. The risk-free rate is 3%, the continuous dividend yield is 1%, and local volatility is $\sigma_{\mathrm{loc}}(S,\tau)=0.2[1+0.2\tanh((S-100)/100)](1+0.1\tau)$.

The diffusion coefficient depends on price and time, and the convection coefficient is $(r-q)S=0.02S$. The spatial domain is $[0,400]$, with a left boundary of 0 and a right boundary of $400e^{-0.01\tau}-100e^{-0.03\tau}$. Compare grids with 201, 301, and 501 spatial nodes and 101, 151, and 251 time points. Also compare two coefficient precomputation strategies on the same medium grid.

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, and finePrice correspond to the coarse, medium, and fine grids. mediumDifference is the absolute difference between the medium and fine grids. cacheDifference is the absolute difference between the "all" and "timeSlice" precomputation strategies on the medium grid.

The maximum difference from the fine grid is approximately 0.00046813 for the medium grid, compared with 0.00138244 for the coarse grid. Both precomputation strategies produce identical results.

Example 4. Down-and-out European call with no dividends.

Price a continuously monitored down-and-out call with no rebate. The current asset price and strike are both 100, the lower barrier is 80, the risk-free rate is 3%, volatility is 20%, and time to maturity is 2 years, with no dividends. The option becomes worthless if the asset touches or falls below 80. If the barrier is never reached, the payoff at maturity is $\max(S-K,0)$.

Restrict the spatial domain to $[80,400]$. Apply a zero Dirichlet condition at the barrier and $400-100e^{-0.03\tau}$ at the upper boundary. Calculate the value at the current asset price of 100 and compare it with the analytical down-and-out call price.

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
)

Example 5. Cash-or-nothing digital call option.

Price a digital call that pays fixed cash at maturity. The current asset price and strike are both 100, the risk-free rate is 3%, volatility is 20%, and time to maturity is 1 year, with no dividends. It pays 5 if the terminal asset price is strictly greater than 100, and 0 otherwise.

The spatial domain is $[0,400]$, with a left boundary of 0 and a right boundary of $5e^{-0.03\tau}$. Use the Crank-Nicolson scheme with rannacherSteps set to 1 to reduce oscillations caused by the discontinuous payoff. Compare the result with the analytical cash-or-nothing call price.

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
)

The relative error on this grid is approximately 0.5742%. Since the payoff is discontinuous at the strike, assess spatial discretization error through grid refinement.

Example 6. Single-barrier call options with a continuous dividend yield.

Price continuously monitored down-and-out and up-and-out calls with no rebate. In both cases, the current asset price and strike are 100, the risk-free rate is 3%, the continuous dividend yield is 1%, volatility is 20%, and time to maturity is 1 year. The convection coefficient is $(r-q)S$.

(1) Down-and-out call: the option becomes worthless when the asset touches or falls below the lower barrier of 80. The spatial domain is $[80,400]$, with a left boundary of 0 and a right boundary of $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) Up-and-out call: the option becomes worthless when the asset touches or exceeds the upper barrier of 130. The spatial domain is $[0,130]$, with zero Dirichlet conditions at both ends. The initial value at the upper barrier is also set to 0 to match the boundary condition. Both cases are compared with their respective analytical prices.

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
)

Example 7. One-touch option with payment at maturity.

Price an upward one-touch option that pays at maturity. The current asset price is 100, the upper barrier is 125, the risk-free rate is 3%, volatility is 20%, and time to maturity is 1 year, with no dividends. If the asset touches 125 during the option's life, it pays 5 at maturity; otherwise, it pays 0.

The spatial domain is $[0,125]$. The terminal payoff is 0 below the barrier and 5 at the barrier. The left boundary is 0, and the barrier boundary is $5e^{-0.03\tau}$, the present value of the now-certain payment at maturity after a touch. Compare the result with the analytical price.

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
)

Example 8. Continuously averaged arithmetic Asian call.

Price a fixed-strike Asian call whose averaging period starts at the valuation time. The current asset price and strike are both 100, the risk-free rate is 5%, volatility is 20%, and both the averaging period and time to maturity are 1 year, with no dividends. The terminal payoff is $\max(A_T/T-K,0)$, where $A_t=\int_0^t S_s\,ds$ is the accumulated price integral and $A_0=0$.

Apply $x=(A-KT)/S$ and $u(S,A,\tau)=(S/T)v(x,\tau)$ to obtain the one-dimensional PDE $v_\tau=\tfrac12\sigma^2x^2v_{xx}+(1-rx)v_x$. Solve on $[-5,5]$ with initial value $\max(x,0)$, a zero left boundary, and the right Neumann condition $v_x=e^{-r\tau}$. Interpolate at $x_0=-1$ and multiply by $S_0/T$ to recover the option price.

For an independent numerical comparison, simulate 200000 Monte Carlo paths with 2000 time steps each and accumulate the price integral using the trapezoidal rule. Set a fixed random seed for reproducibility and report the Monte Carlo standard error.

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 is the Monte Carlo estimate, mcStdError is its sampling standard error, and absoluteDifference is its absolute difference from the PDE price. The difference is approximately 1.46 standard errors.

Example 9. European call under the Heston stochastic volatility model.

Price a European call driven by the asset price and its instantaneous variance. The current asset price and strike are both 100, the risk-free rate is 5%, and time to maturity is 1 year, with no dividends. The initial variance is $v_0=0.04$, long-run variance is $\theta=0.04$, mean-reversion speed is $\kappa=2$, volatility of variance is $\xi=0.3$, and price-variance correlation is $\rho=-0.7$.

The code uses $x=\ln S$ and variance $v$ as the two coordinates, with prices in $[1,2000]$ and variances in $[0,1]$. The terminal payoff is $\max(e^x-100,0)$. The lower and upper price boundaries are 0 and $2000-100e^{-0.05\tau}$, respectively. Zero Neumann conditions approximate both variance boundaries.

Use the MCS ADI scheme for the two-dimensional PDE with a mixed derivative and interpolate at [log(100), 0.04]. Compare the result with the reference price of 10.394218565150163 obtained by numerical integration of the Heston characteristic function for the same parameters.

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
)

Example 10. Two-asset basket call option.

Price a European call on an equally weighted basket of two correlated assets. Both current asset prices are 100, both weights are 0.5, the strike is 100, the risk-free rate is 5%, and time to maturity is 1 year, with no dividends. The volatilities are 20% and 25%, correlation is 0.3, and the terminal payoff is $\max(0.5S_1+0.5S_2-100,0)$.

Both asset price domains are $[0,400]$. When one asset price is zero, the boundary value is the Black-Scholes call price for the weighted value of the other asset. At both upper price boundaries, impose Neumann conditions with an outward normal derivative of 0.5. The xy coefficient represents the mixed derivative arising from asset correlation. Use the MCS ADI scheme to solve the PDE.

Interpolate at asset prices of 100 and 100, and compare the result with a Monte Carlo price from 1000000 pairs of correlated terminal prices. Set a fixed random seed and report the sampling standard error.

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 is the absolute difference between the PDE price and the Monte Carlo estimate, approximately 0.31 Monte Carlo standard errors in this example.

Example 11. Solve the 1D heat equation $u_t=0.2u_{xx}$ on $[0,1]$, with zero Dirichlet conditions at both endpoints and initial condition $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

Example 12. Solve the 2D heat equation $u_t=0.05(u_{xx}+u_{yy})$ on $[0,1]\times[0,1]$, with zero Dirichlet conditions on all four boundaries.

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

Example 13. Solve a two-dimensional Poisson-type problem. The boundary vector is ordered as xMin, xMax, yMin, yMax: xMin and yMin use Robin conditions, while xMax and yMax use Neumann conditions.

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

Example 14. On a DolphinDB Server with JIT support, set the jit field in method to true to request JIT compilation for coefficient, initial-condition, and boundary callbacks. Define named user-defined functions, then pass their names to the corresponding arguments.

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

Example 15. Set time-dependent Dirichlet boundaries for $\frac{\partial u}{\partial t}=0.2\frac{\partial^2u}{\partial x^2}+1$ on [0,1], with initial condition $u(x,0)=x$ and boundaries $u(0,t)=t$ and $u(1,t)=1+t$. Pass a def(t) function directly for the left endpoint and use a dictionary for the right endpoint to illustrate two equivalent representations. The initial values match the boundaries at t=0. The exact solution is $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]

Example 16. Set nonzero Neumann boundaries and combine Dirichlet and Neumann boundaries. Solve $\frac{\partial u}{\partial t}=0.2\frac{\partial^2u}{\partial x^2}-0.4$ on [0,1.5] with initial condition $u(x,0)=x^2+3x+5$. The exact solution is stationary. The Neumann value specifies the outward normal derivative: $-\left.\frac{\partial u}{\partial x}\right|_{x=0}=-3$ at the left endpoint and $\left.\frac{\partial u}{\partial x}\right|_{x=1.5}=6$ at the right endpoint. This example uses a nonuniform grid; the default second-order boundary differences require at least 4 nodes. The second solve fixes the left endpoint at 5 and retains the right Neumann boundary.

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]

Example 17. Use functions for alpha, beta, and value in Robin boundary dictionaries. Solve $\frac{\partial u}{\partial t}=0.2\frac{\partial^2u}{\partial x^2}+1$ with initial condition $u(x,0)=x$ and exact solution $u(x,t)=x+t$. At both endpoints, set $\alpha(t)=1+t$ and $\beta(t)=0.5+0.1t$. From $\alpha(t)u+\beta(t)\frac{\partial u}{\partial n}=\mathrm{value}(t)$, the left endpoint requires $\mathrm{value}(t)=(1+t)t-(0.5+0.1t)$, and the right endpoint requires $\mathrm{value}(t)=(1+t)^2+(0.5+0.1t)$. The minus and plus signs follow the outward normal directions at the left and right endpoints.

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]

Related functions: pdeCoeff1D, pdeCoeff2D, pdeDomain1D, pdeDomain2D, pdeGrid1D, pdeInfo