# How to implement FEM for a 2D PDE with variable coefficients

Consider the following 2D linear time-dependent PDE with variable coefficients: $$$$\label{Bates_equation} \begin{split} \frac{\partial u(s,v,\tau)}{\partial\tau} &=\frac{1}{2}v s^{2}\frac{\partial^{2}u(s,v,\tau)}{\partial s^{2}}+\frac{1}{2}\sigma^{2}v\frac{\partial^{2}u(s,v,\tau)}{\partial v^{2}}\\ &\quad +\rho\sigma v s\frac{\partial^{2}u(s,v,\tau)}{\partial s\partial v}\\ &\quad +(r-q-\lambda\xi)s\frac{\partial u(s,v,\tau)}{\partial s}+\kappa(\theta-v)\frac{\partial u(s,v,\tau)}{\partial v}\\ &\quad-(r+\lambda)u(s,v,\tau), \end{split}$$$$ where $$0\leq s \leq s_{max}=3E$$, $$0\leq v\leq 1$$, E=100 and some other suitable values. To solve this PDE using the FEM, I employ the answer given by user21 at FEM differentiation matrices, to construct the FEM differentiation matrices and then obtain the corresponding system of ODEs that must be solved to find the final solution. I did this as follows:

ClearAll["Global*"];
Needs["NDSolveFEM"]
FiniteElementDerivative[order : {__Integer}, mesh_ElementMesh] /;
1 <= Length[order] <= 3 :=
Block[{dim, nr, vd, sd, mdata, ccoef, pos, dcoef, cdata},
dim = Length[order];
nr = ToNumericalRegion[mesh];
vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u},
Table[Unique[X], {dim}]}];
sd = NDSolveSolutionData[{"Space"} -> {nr}];
mdata = InitializePDEMethodData[vd, sd];
ccoef = ConstantArray[0, dim];
pos = Flatten[Position[order, 1]];
ccoef[[pos]] = 1;
dcoef = ConstantArray[0, dim];
pos = Flatten[Position[order, 2]];
dcoef[[pos]] = 1;
dcoef = DiagonalMatrix[dcoef];
(*"Pure ConvectionCoefficients" will trigger a warning*)
Quiet[cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{dcoef}},
"ConvectionCoefficients" -> {{ccoef}}],
{InitializePDECoefficients::femcscd}];
DiscretizePDE[cdata, mdata, sd]]


And then write:

m = 32; n = 32; size = m*n;
TT = 1.; r = 0.025;
e = 100.; q = 0; \[Sigma] = 0.3; \[Kappa] = 1.5; \[Theta] = 0.04; \
\[Lambda] = 0; \[Zeta] = 0; \[Rho] = -0.9;
\[Rho]1 = Sqrt[1 - \[Rho]^2]; \[Omega] = r - q - \[Lambda]*\[Zeta];
xsmin = 0; xsmax = 3 e; ysmin = 0; ysmax = 1.0;

hx = (xsmax - xsmin)/(m - 1); nx = N@Range[xsmin, xsmax, hx]; nx1 =
Partition[nx, 1];
hy = (ysmax - ysmin)/(n - 1); ny = N@Range[ysmin, ysmax, hy];
origrid = Flatten[Outer[List, nx, ny], 1];

Idx = SparseArray[{{i_, i_} -> 1.}, {m, m},
0]; xx = (SparseArray@DiagonalMatrix@nx);
Idy = SparseArray[{{i_, i_} -> 1.}, {n, n},
0]; yy = (SparseArray@DiagonalMatrix@ny);
XX = KroneckerProduct[xx, Idy]; Y = KroneckerProduct[Idx, yy];
Id = KroneckerProduct[Idx, Idy];

mesh = ToElementMesh[origrid];
mesh["Wireframe"];

{dXmatFEM, d2XmatFEM} =
FiniteElementDerivative[#, mesh]["StiffnessMatrix"] & /@ {{1,
0}, {2, 0}};
Map[MatrixPlot, {dXmatFEM, d2XmatFEM}];
{dYmatFEM, d2YmatFEM} =
FiniteElementDerivative[#, mesh]["StiffnessMatrix"] & /@ {{0,
1}, {0, 2}};
dXYmatFEM = (FiniteElementDerivative[#, mesh][
"StiffnessMatrix"] & /@ {{1, 1}})[[1]];
Map[MatrixPlot, {dYmatFEM, d2YmatFEM, dXYmatFEM}];

mat = SparseArray@Simplify[(1/2 Y.XX^2).d2XmatFEM
+ ((\[Rho]*\[Sigma])*(XX.Y)).(dXYmatFEM) + (1/
2 \[Sigma]^2 Y).d2YmatFEM
+ (\[Omega]*
XX).(dXmatFEM) + (\[Kappa]*(\[Theta]*Id -
Y)).(dYmatFEM) - (r + \[Lambda]) Id];
U[t_] = Flatten@Table[Subscript[u, i, j][t], {i, 1, m}, {j, 1, n}];
(*Initial condition*)
payoff = Flatten@Table[Max[nx[[i]] - e, 0], {i, 1, m}, {j, 1, n}];
eqns = Thread[D[U[t], t] == mat.U[t]];

(*Boundaries*)
sf = 0; bc1 = Table[Subscript[u, 1, j][t] == 0., {j, 1, n}];
Table[bc11[l] = Map[D[#, t] + sf # &, bc1[[l]]], {l, 1, Length[bc1]}];
sf = 0; bc2 =
Table[Subscript[u, m, j][t] == xsmax*Exp[-q t] - e*Exp[-r t], {j, 1,
n}];
Table[bc22[l] = Map[D[#, t] + sf # &, bc2[[l]]], {l, 1, Length[bc2]}];

sf = 0; bc4 =
Table[Last[
NDSolveFiniteDifferenceDerivative[1, ny,
Take[U[t], {(k - 1)*n + 1, (k) n}]],
"DifferenceOrder" -> 2] == 0, {k, 2, m - 1}];
Chop@Table[
bc44[l] = Map[D[#, t] + sf # &, bc4[[l]]], {l, 1, Length[bc4]}];

Table[eqns[[i]] = bc11[i], {i, 1, Length[bc1]}];
Table[eqns[[size - i + 1]] = bc22[Length[bc2] - i + 1], {i,
Length[bc2], 1, -1}];
Table[eqns[[(i + 1) n]] = bc44[i], {i, 1, Length[bc4]}];

For[k = 2, k <= m - 1, k++,
bo1 = Drop[Take[U[t], {(k - 1)*n + 1, (k) n}], -1];
eqns[[(k*n)]] = (Chop@Simplify[
Last@Table[eqns[[(k*n)]] = eqns[[(k*n)]] /. (D[bo1[[i]], t] ->
Last@eqns[[((k - 1)*n + 1) + (i - 1)]]), {i, 1,
Length@bo1}]
]);
coef1 = Normal@CoefficientArrays[eqns[[k*n]], U[t]];
coef2 = Coefficient[First[coef1], D[U[t][[k*n]], t]];
eqns[[k*n]] =
D[U[t][[(k*n)]], t] == Chop@(-(1/coef2) Last@coef1).U[t];
];

vec0 = SparseArray[{i_} -> 0, size];
mat01 = Table[vec0, {i, 1, n}];
mat02 = Table[-Last@CoefficientArrays[eqns[[i]], U[t]], {i,
n + 1, (m - 1) n}];
mat03 = ArrayFlatten[{{mat01}, {mat02}, {mat01}}];
mat03 // MatrixPlot;
vec1 = Last@eqns[[m*n]];
vec2 = SparseArray[{i_} -> 0, ((m - 1)*n)];
vec3 = SparseArray[{i_} -> vec1, n];
vec4[t_?NumericQ] = N@Join[vec2, vec3];
mat03 = SparseArray[Chop@mat03];

(*SOLVING SYSTEM OF ODES*)
Monitor[lines =
NDSolve[{D[v[t], t] == mat03.v[t] + vec4[t],
v[0] == initc[[All, 2]]}, v[t], {t, 0, TT},
AccuracyGoal -> 5, PrecisionGoal -> 5,
Method -> {"FixedStep", "StepSize" -> .0001,
Method -> "ExplicitEuler"},
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
monitor]; // AbsoluteTiming

s = v[t] /. lines[[1]]; tss = s[[0]][[3]][[1]];
Print["The number of total time steps = ", Length@tss];

s /. {t -> 0};
set0 = Table[Flatten@{origrid[[i]], %[[i]]}, {i, 1, size}];
ListPlot3D[set0, AxesLabel -> {"s", "v", "u"}, PlotRange -> All,
ImageSize -> 400]

sol0 = s[[0]][[4]]; sol1 = s[[0]][[4]][[Length@sol0]];
sol2 = sol1[[1]];
list1 = Table[Flatten@{origrid[[i]], sol2[[i]]}, {i, 1, size}];

T12 = Map[Last, list1];
set1 = Table[Flatten@{origrid[[i]], T12[[i]]}, {i, 1, size}];
ListPlot3D[set1, AxesLabel -> {"s", "v", "u"}, PlotRange -> All,
ImageSize -> 400]

v0 = 0.04;
f = Interpolation@list1;
Print["The error is = ", ScientificForm@Abs[8.894869 - f[100, v0]]];


Unfortunately, this does not converge to the solution as the last line showing by not approaching zero when I increase the number of nodes, $$m$$ and $$n$$.

Can anyone give some hints how I can impose the variable coefficients in the final system matrix? Maybe I am making a mistake at that part! Or the procedure for constructing the matrix in the system of discretized ODEs is different in FEM methodology?

The boundary conditions are: $$$$\label{CH3Eq4} \begin{split} &u(s,v,\tau)\simeq 0,\qquad s\rightarrow0,\\ &u(s,v,\tau)\simeq s_{\text{max}}\exp{(-q\tau)}-E\exp{(-r\tau)},\qquad s\rightarrow+\infty,\\ &\frac{\partial u(s,v,\tau)}{\partial v}\simeq 0,\qquad v\rightarrow+\infty. \end{split}$$$$

A sample numerical solution for PDE must look like the following:

• How you determined the the error is of 8.894869 - f[100, v0]? Feb 11 '21 at 22:02
• We have the reference solution from a published work which says that at a specific point (100, 0.04) of the domain, the solution is 8.8948. Accordingly, after I compute the numerical solution by solving the system of linear ODEs, and save it in $f$, I can find the final error. The point here is, if we use the FD method and fill up the differentiation matrices with FD weights, the last part of the code would be similar and everything goes fine and you can see the convergence. But for the FEM, I do not know how to implement the FEM differentiation matrices to construct the final system of ODEs! Feb 12 '21 at 8:28
• Could you give a link to the paper you mentioned above? Feb 12 '21 at 10:15
• Please see Table 3 (Case 1) in the PDF that can be downloaded in the following link: semanticscholar.org/paper/… Feb 12 '21 at 13:03
• Boundary conditions at infinity look useless. We can try to solve problem without boundary conditions, only with initial data as Monitor[lines = NDSolve[{eqns, initc}, U[t], {t, 0, TT}, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor];. Then we can try to implement bc at v=0. Feb 13 '21 at 22:03

We can solve this problem using FEM and implicit Euler in a Case 1, 2, 4 from the report Option Pricing under a Heston Volatility model using ADI schemes by Jieshun Luo, Qi Wang, Nestor Carbayo. Case 3 looks very tricky, so I can't handle it with my code. First we map region $$0\le s\le 800, 0\le v\le 5$$ to unit square by substitution $$x=s/L_x,y=v/L_y$$ with $$L_x=800, L_y=5$$ and define mesh as follows (in this code we use notation from above)

Needs["NDSolveFEM"]

TT = 1.; r = 0.025; eps = 10^-3;
e = 100.; q = 0; \[Sigma] = 0.3; \[Kappa] = 1.5; \[Theta] = 0.04; \
\[Lambda] = 0; \[Zeta] = 0; \[Rho] = -0.9; Lx = 8 e; Ly = 5; xc =
e/Lx; yc = .04/Ly; p0 = 8.894869;
\[Rho]1 = Sqrt[1 - \[Rho]^2]; \[Omega] = r - q - \[Lambda]*\[Zeta];
reg = Rectangle[{eps, eps}, {1, 1}]; f =
Function[{vertices, area},
Block[{x, y}, {x, y} = Mean[vertices];
If[Abs[x - xc] <= xc/4 || (y <= xc/4 && 0 < x < 1), area > 0.00005,
area > 0.0005]]]; mesh =
ToElementMesh[reg, AccuracyGoal -> 5, PrecisionGoal -> 5,
MeshRefinementFunction -> f];mesh1 = ToElementMesh[Line[{{eps}, {1}}], MaxCellMeasure -> .001];
mesh["Wireframe"]


Second step, we define equation, initial and boundary conditions in variables x,y as

eq = (u[x, y] - U[i - 1][x, y])/
dt - (Ly/2 x^2 y D[u[x, y], x, x] + ((\[Rho]*\[Sigma])*x y) D[
u[x, y], x, y] + (1/2 \[Sigma]^2 y/Ly) D[u[x, y], y,
y] + (\[Omega]*x) D[u[x, y],
x] + (\[Kappa]*(\[Theta] - y Ly)/Ly) D[u[x, y],
y] - (r + \[Lambda]) u[x, y]);
U[0][x_, y_] := Lx If[x - xc > 0, x - xc, 0];
Ub[0] = U[0][x,
eps]; bc = {DirichletCondition[u[x, y] == 0, x == eps],
DirichletCondition[u[x, y] == Lx x Exp[-q i dt], y == 1],
DirichletCondition[u[x, y] == Ub[i - 1], y == eps]};
bc1 = NeumannValue[y Ly Lx/2 Exp[- q i dt], x == 1];


Note, that we need to prepare code for implicit Euler and there is also additional equation describing unknown function Ub (last equation in section 2.2 of the report), finally we have

nt = 100; dt = 1/nt;
Do[U[i] =
NDSolveValue[{eq == bc1, bc}, u, Element[{x, y}, mesh]] // Quiet;
sol = NDSolveValue[{-D[w[t, x], t] + (\[Omega]*x) D[w[t, x],
x] + (\[Kappa]*(\[Theta])/Ly) Derivative[0, 1][U[i]][x,
eps] - (r + \[Lambda]) w[t, x] == 0,
DirichletCondition[w[t, x] == 0, x == eps],
w[(i - 1) dt, x] == Ub[i - 1],
DirichletCondition[w[t, x] == U[i][1, eps], x == 1]},
w[t, x], {t, (i - 1) dt, i dt}, Element[{x}, mesh1]] // Quiet;
Ub[i] = sol /. {t -> i dt}, {i, 1, nt}] // AbsoluteTiming


It takes about 1 min to get final step, and solution on this step at t=1 looks like

Plot3D[U[nt][x/Lx, y/Ly], {x, eps Lx, Lx}, {y, eps Ly, Ly - .2},
PlotRange -> All, AxesLabel -> {"s", "v"}, Mesh -> None,
ColorFunction -> "Rainbow"]


In the reference point s=100, v=0.04 we have

U[nt][xc, yc]

Out[]= 8.86571


We can compare it with p0 = 8.894869 to estimate error. Note, the boundary condition at v=0 has a form $$u(t,s,0)=w(t, s)$$ and $$\frac {\partial w}{\partial t}=\omega x \frac {\partial w}{\partial s}+\kappa\theta \frac {\partial u}{\partial v}|_{v=0}-(r + \lambda) w$$ Therefore we need to solve this PDE to use w as a boundary condition for u. It is why we need to implement implicit algorithm for u.

Update 1. Code for Case 2 from report. This is different from above since we have used notation from report.

Needs["NDSolveFEM"]
TT = 1.; rd = 0.01; rf = 0.04; eps = 10^-3;
e = 100.; \[Sigma] = 0.04; \[Kappa] = 3; \[Rho] = 0.6; \[Eta] = 0.12; \
Lx = 8 e; Ly = 5; xc = e/Lx; yc=.04/Ly;
reg = Rectangle[{eps, eps}, {1, 1}]; f =
Function[{vertices, area},
Block[{x, y}, {x, y} = Mean[vertices];
If[Abs[x - xc] <= xc/4 || (y <= xc/4 && 0 < x < 1) ||
x >= 1 - xc/4 || y >= 1 - xc/4, area > 0.000075,
area > 0.00075]]]; mesh =
ToElementMesh[reg, MeshRefinementFunction -> f];

mesh1 = ToElementMesh[Line[{{eps}, {1}}], MaxCellMeasure -> .00075];
eq = (u[x, y] - U[i - 1][x, y])/
dt - (Ly/2 x^2 y D[u[x, y], x, x] + ((\[Rho]*\[Sigma])*x y) D[
u[x, y], x, y] + (1/2 \[Sigma]^2 y/Ly) D[u[x, y], y,
y] + (rd - rf) x D[u[x, y],
x] + (\[Kappa]*(\[Eta] - y Ly)/Ly) D[u[x, y], y] - rd u[x, y]);
U[0][x_, y_] := Lx If[x - xc > 0, x - xc, 0];
Ub[0] = U[0][x,
eps]; bc = {DirichletCondition[u[x, y] == 0, x == eps],
DirichletCondition[u[x, y] == Lx x Exp[-rf i dt], y == 1],
DirichletCondition[u[x, y] == Ub[i - 1], y == eps]};

nt = 100; dt = 1/nt;
Do[U[i] =
NDSolveValue[{eq == NeumannValue[Lx Ly/2 y Exp[- rf i dt], x == 1],
bc}, u, Element[{x, y}, mesh]] // Quiet;
sol = NDSolveValue[{-D[w[t, x], t] + (rd - rf) D[w[t, x],
x] + (\[Kappa]*\[Eta]/Ly - eps) Derivative[0, 1][U[i]][x,
eps] - rd w[t, x] == 0,
DirichletCondition[w[t, x] == 0, x == eps],
w[(i - 1) dt, x] == Ub[i - 1]}, w[t, x], {t, (i - 1) dt, i dt},
Element[{x}, mesh1]] // Quiet;
Ub[i] = sol /. {t -> i dt}, {i, 1, nt}] // AbsoluteTiming


Visualization of mesh and solution

{mesh["Wireframe"],
Plot3D[U[nt][x/Lx, y/Ly], {x, eps Lx, Lx}, {y, eps Ly, Ly - .2},
PlotRange -> All, AxesLabel -> {"s", "v"}, Mesh -> None,
ColorFunction -> "Rainbow"]}


Finally we compute solution in the reference point and estimate error

U[nt][xc, yc]

Out[]= 10.5396

(10.541784 - U[nt][xc, yc])/10.541784

Out[]= 0.000210381

• Thank you. Your implementation is useful, however, the only thing is that your final result is a bit different in terms of accuracy because of boundary conditions. I think the boundary at one side of the domain is not imposed well. Maybe it is inconsistent with the methodology of FEM and yields such figures and not accurate final result. Can you please take a look at the boundary conditions (y tends to zero and one) and also see the sample solution picture that I now edited and included above? Or maybe it would be better if you do the transformation back and then plot the final solution. Feb 14 '21 at 21:06
• Do you understand that there are Case 1 and Case 2 solved (last code is Case 2 solution)? I have used equations (5)-(8) from report as initial and boundary conditions. And also last equation from section 2.2. My final solution has error of $10^{-3}$ in a reference point. Your boundary condition and region $0\le s \le 3 e, 0\le v \le 1$ not compatible with data from report (they used region $0\le s \le 8 e, 0\le v \le 5$). Feb 14 '21 at 21:34
• Thanks. $s_{\max}$ and $v_{\max}$ can be the values that I wrote or the values in the report that you used. However, OK. Can you please draw the final solution curve (at least for Case 1) in the original domain instead of the transformed domain? Feb 14 '21 at 21:58
• And what's the reason of using an implicit ODE solver? What do you think about using other solvers? Feb 14 '21 at 21:59
• Actually in this code the implicit PDE solver used. We solve u[x,y] on every step on time. We can't implement bc at y=0 directly in NDSolve[] since it can't accept bc in a form of PDE for w[t,x] as in my code. Feb 14 '21 at 23:01