# Generate coefficient array from general formula of linear equation system

This is a problem coming out in the implementation of finite difference method (FDM).

Here is a toy example. Suppose we want to solve the boundary value problem (BVP)

$$y''(x)=\sin(x),\ y(0)=0,\ y(1)=0$$

eq = y''[x] == Sin[x];
bc = {y[0] == 0, y[1] == 0};


with FDM, then we can program as follows. First, write down the corresponding difference formula:

Clear@dx;
generalformula = (y[x + dx] - 2 y[x] + y[x - dx])/dx^2 == Sin[x];


Then define the grid size and generate difference equations at every grid point:

points = 25; domain = {0, 1};
dx = Subtract @@ domain/(1 - points);
var = Table[y[x], {x, ##, dx}] & @@ domain;
differenceeq = Table[generalformula, {x, domain[[1]] + dx, domain[[-1]] - dx, dx}]~Join~bc;


Finally solve the difference equations with e.g. NSolve:

solrule = NSolve[differenceeq, var][[1]]; // AbsoluteTiming
(* {0.00462427, Null} *)


BTW, a comparision with the analytic solution:

asol = DSolveValue[{eq, bc}, y, x]
Plot[asol[x], {x, ##}, PlotRange -> All]~Show~
ListPlot[solrule[[All, -1]], DataRange -> {##}] & @@ domain


Alternatively, we can extract the coefficient array from differenceeq and solve it with LinearSolve, as illustrated in the document of CoefficientArrays. This approach is faster than the one using NSolve:

{barray, marray} = CoefficientArrays[differenceeq, var]; // AbsoluteTiming
(* {0.000614324, Null} *)
sol = LinearSolve[marray // N, -barray]; // AbsoluteTiming
(* {0.000488251, Null} *)


However, this isn't the end. As already mentioned, we've generated difference equation at every grid point and then extract the coefficient array from the symbolic system in the code above. This procedure can be rather memory-demanding and slow when points is large. Is it possible to generate barray and marray directly from generalformula and bc programatically, in a memory-conserving and not-that-slow way?

If the toy example above is too simple to illustrate the underlying problem, the following is a more complicated one attempting to solve the problem mentioned here:

r = .8; glass = (1.45)^2; air = 1.; k = (2 π)/1.55;
n2 = Function[{x, y}, #] &@If[x^2 + y^2 <= r^2, glass, air];
ClearAll[fw, bw, delta]
SetAttributes[#, HoldAll] & /@ {fw, bw};

fw@D[expr_, x_] :=
With[{Δ = delta@ToString@x},
Subtract @@ (expr /. {{x -> x + Δ}, {x -> x}})/Δ]
bw@D[expr_, x_] :=
With[{Δ = delta@ToString@x},
Subtract @@ (expr /. {{x -> x}, {x -> x - Δ}})/Δ]

delta[a_ + b_] := delta@a + delta@b
delta[k_. delta[_]] := 0
grad = Function[expr, fw@D[expr, #] & /@ {x, y}];
div = Function[expr, Total@MapThread[bw@D@## &, {expr, {x, y}}]];
elst = e[#][x, y] & /@ Range[2];
discretelhs =
With[{n2 = n2[x, y]},
div@grad@elst + (n2 k^2) elst - grad[div@elst - 1/n2 div[n2 elst]]] // FullSimplify;
points = 100; L = 2; domain = {-L, L}; grid = Array[# &, points, domain];
del = #[[2 ;; -2]] &;
delta["x"] = delta["y"] = -Subtract @@ domain/(points - 1);
vare = Outer[e[#][#2, #3] &, Range@2, del@grid, del@grid, 1] // Flatten;
discretelhslst =
Flatten@Transpose[
Table[discretelhs, {x, del@grid}, {y, del@grid}], {2, 3, 1}]; // AbsoluteTiming
{barray2, marray2} = CoefficientArrays[discretelhslst, vare]; // AbsoluteTiming
{val, vec} =
Eigensystem[marray2, -6,
Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
MapThread[ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All,
ColorFunction -> "AvocadoColors"] &, {mat, val}]


ByteCount[marray2]/1024^2. "MB"
(* 1.60567 "MB" *)

ByteCount[discretelhslst]/1024^2. "MB"
(* 44.9877 "MB" *)


The corresponding problem is, how to arrive at marray2 from discretelhs without generating discretelhslst, in an efficient enough way?

• Would a pre-processed form that is only really applicable for FDM-type matrices be acceptable? I can think up a good way to do that. Commented Jun 19, 2019 at 7:36
• @b3m2a1 I'm only interested FDM here. What do you mean by "pre-processed form"? If it means something like standard form, I think a requirement is, it should not be too hard to transform a difference formula to that standard form. Commented Jun 19, 2019 at 7:46
• I mean that I take what you supply as this FDM expression and extract the weights from it. Then I can make the FDM matrix easily using sparse methods. Commented Jun 19, 2019 at 7:47
• @b3m2a1 Sounds like a proper answer for my question. Commented Jun 19, 2019 at 7:49
• @Alex Again, it's surprising that removing //N improves the performance. Actually I add it to handle grid because in certain cases complicated symbolic grid can slow down the code, but I didn't expect inexact grid is more complicated in this case. Edited, thx for pointing out. However, Block should not be replaced by Module, or the Dirichlet b.c. won't be set. Just try expr = f[4]; Module[{f}, f[x_] =x^2; expr]. The code still produces reasonable result because b.c. isn't that important for this problem. Another example is here: mathematica.stackexchange.com/q/163486/1871 Commented Jun 29, 2019 at 2:52

### Partial ND

It's possible to again extract all this weight data and stuff in the higher dimensional cases using effectively the same procedure as for the 1D:

fdmWeightDataND~SetAttributes~HoldAll;
fdmWeightDataND[a_, dep_Symbol, indeps : {__Symbol}, hs : {__Symbol}] :=

Block[{dep}, Block[indeps, Block[hs,
Module[{weights, order, denomWeight, denom, num, powed, pows, minpow,
pts, x}, denom = Denominator[a];
(*assumes the mesh spacing is in there at different powers*)

order = Exponent[denom, #] & /@ hs;
(*pull the power on the mesh spacing*)

denomWeight = If[order > 0, Coefficient[denom, hs^order], denom];
num = Numerator[a];
{powed, pows} =
Reap[
num /.
dep[args___]?(Not@*FreeQ[Alternatives @@ indeps]) :>
dep @@
Sow[
Coefficient[#, #2] &,
{
{args},
hs
}
]
]
];
minpow = Min[pows];
weights =
pows[[1]],
Map[Coefficient[powed /. (dep @@ #) -> x, x] &, pows[[1]]]
];
<|
"Weights" -> weights,
"Orders" -> order
|>
]
]
]
];


This is less automate-able, though, since you need to specify what your f, independent variables, and mesh spacings are. Also without knowing how you'd like to flatten your nD tensor of grid points it's impossible to automatically build the FD matrix. But at a minimum I can give you weights and indices for the weights:

(*didn't want to Wikipedia the nine-point stencil so most of this is fake*)

simpleFormula = (
f[x - xh, y] + f[x + xh, y] +
f[x, y - yh] + f[x, y + yh] - 4 f[x, y] +
1/2 f[x + xh, y - yh] + 1/2 f[x - xh, y - yh] +
1/2 f[x + xh, y + yh] + 1/2 f[x - xh, y + yh]
)/(xh^2*yh^2);

weights2D = fdmWeightDataND[simpleFormula, f, {x, y}, {xh, yh}]

<|"Weights" -> <|{0, 0} -> -4, {0, -1} -> 1, {0, 1} -> 1, {-1, 0} ->
1, {-1, -1} -> 1/2, {-1, 1} -> 1/2, {1, 0} -> 1, {1, -1} -> 1/2, {1, 1} ->
1/2|>, "Orders" -> {2, 2}|>


Basically now you'd pick how to order your grid points and build a matrix from that. Here's one way to do this based on having the second index be the slowest moving (note that I'm only like 50% confident I built this correctly):

FDM2D[weightData_, {hx_, hy_}, {nptsx_, nptsy_}] :=

Module[{xps, yps, weights = weightData["Weights"],
orders = weightData["Orders"]},
xps = Length@DeleteDuplicates@Keys[weights][[All, 1]];
yps = Length@Sort@DeleteDuplicates@ Keys[weights][[All, 2]];
SparseArray[
Flatten@
KeyValueMap[
With[{ix = #[[1]], iy = #[[2]], v = #2},
Map[
Band[
If[ix >= 0, {1, 1 + ix}, {1 - ix, 1}] + {#, #} +
If[iy >= 0, {0, iy*nptsx}, {-iy*nptsx, 0}],
{nptsx*nptsy, nptsx*nptsx},
{nptsx, nptsx}
] -> v &,
Range[nptsx - Abs[ix]] - 1
]
] &,
weights
],
{nptsx*nptsy, nptsx*nptsx}
]/(Times @@ ({hx, hy}^orders))
]


Here's what that outputs:

FDM2D[weights2D, {1, 1}, {5, 5}] // MatrixPlot


In the case of no coupling across degrees of freedom this looks correct I believe... (but can't be sure since I usually do my Laplacian via a Kronecker sum):

laplacian = (
f[x - xh, y] + f[x + xh, y] +
f[x, y - yh] + f[x, y + yh] - 4 f[x, y]
)/(xh^2*yh^2);

FDM2D[fdmWeightDataND[laplacian, f, {x, y}, {xh, yh}], {1, 1}, {10,
10}] // MatrixPlot


### Full Thing for 1D

Here's a method that automates this extraction for 1D FDM stuff:

fdmWeightData~SetAttributes~HoldAll;
fdmWeightData[Equal[a_, b_], dep : {__Symbol}, indep : {__Symbol}] :=
fdmWeightData[a, dep, indep];
fdmWeightData[a_, dep_Symbol, indep_Symbol, h_Symbol] :=

Block[{dep, indep, h},
Module[{weights, order,
denomWeight, denom,
num, powed, pows,
minpow, pts
},
denom = Denominator[a];
(* assumes the mesh spacing is in there *)

order = Exponent[denom, h];
(* pull the power on the mesh spacing *)

denomWeight = Coefficient[ denom, h^order];
num = Numerator[a];
{powed, pows} =
Reap[
num /.
dep[expr_?(Not@*FreeQ[indep])] :>
dep[Sow@Coefficient[expr, h]]
];
minpow = Min[pows];
powed = powed /. dep[n_] :> dep^(n - minpow);
weights = (CoefficientList[powed, dep])/denomWeight;
pts = Exponent[powed, dep, List] + minpow;
<|
"Weights" -> weights,
"Points" -> pts,
"Order" -> order
|>
]
];
fdmWeightData[Equal[expr_, b_]] :=
fdmWeightData[expr];
fdmWeightData[expr_] :=
Module[
{
vars = Variables[expr],
dep,
indep,
h
},
h = FirstCase[vars, _Symbol];
{dep, indep} =
Replace[
FirstCase[vars, a_Symbol[b_Symbol] :> {a, b}, None],
None :>
a_Symbol[e_] :> FirstCase[e, Except[h, _Symbol], None, Infinity]
] ;
With[
{dep = dep, indep = indep, h = h},
fdmWeightData[expr, dep, indep, h]
]
];

fdmWeightData[(y[x + dx] - 2 y[x] + y[x - dx])/dx^2 == Sin[x]]

<|"Weights" -> {1, -2, 1}, "Points" -> {-1, 0, 1}, "Order" -> 2|>


Here's a helper on that:

fdmMat[weights_List, pts_List, order_List, h_, n_] :=
SparseArray[
If[# >= 0,
Band[{1, # + 1}] -> #2,
Band[{-# + 1, 1}] -> #2
] &,
{pts, weights/(h^order)}
],
{n, n}
];
fdmMat[a_Association, h_, n_] :=

fdmMat[Sequence @@ Lookup[a, {"Weights", "Points", "Order"}], h, n];
fdmMat[a_Association _, n_] :=
fdmMat[a, 1, n];
FDMMatrix[a_Association][n_, h_: 1] :=
fdmMat[a, h, n];
FDMMatrix[expr : Except[_Pattern | _Association]] :=

FDMMatrix[fdmWeightData[expr]]
FDMMatrix[expr_, dep_, indep_, h_] :=

FDMMatrix[fdmWeightData[expr, dep, indep, h]]

fdmat = FDMMatrix[(y[x + dx] - 2 y[x] + y[x - dx])/dx^2]

FDMMatrix[<|"Weights" -> {1, -2, 1}, "Points" -> {-1, 0, 1}, "Order" -> 2|>]


Then we can make the matrix fast:

fdmat[100] // MatrixPlot


Doing the nD extension will take a bit of thought, though.

### Partial Solution to Actual Problem

I'm not sure how your actual problem works, but here's a way we can try to build an FD matrix for it using this stuff. The main thing we'll do is break it down into portions with the same dx and dy terms after doing some symbolic pre-simplification. Here's my new setup (same as yours but entirely symbolic):

Clear[r, glass, air, k]
n2 = (Function[{x, y}, #] &@(SimplifyPWToUnitStep@
PiecewiseExpand@If[x^2 + y^2 <= r^2, glass, air]));
ClearAll[fw, bw, delta]
SetAttributes[#, HoldAll] & /@ {fw, bw};

fw@D[expr_, x_] :=
With[{\[CapitalDelta] = delta@ToString@x},
Subtract @@ (expr /. {{x -> x + \[CapitalDelta]}, {x -> x}})/\[CapitalDelta]]
bw@D[expr_, x_] :=
With[{\[CapitalDelta] = delta@ToString@x},
Subtract @@ (expr /. {{x -> x}, {x -> x - \[CapitalDelta]}})/\[CapitalDelta]]

delta[a_ + b_] := delta@a + delta@b
delta[k_. delta[_]] := 0
grad = Function[expr, fw@D[expr, #] & /@ {x, y}];
div = Function[expr, Total@MapThread[bw@D@## &, {expr, {x, y}}]];
elst = e[#][x, y] & /@ Range[2];
baseDiscreteLHS =
With[{n2 = n2[x, y]},


Now I'll prep this stuff so it's in a form my code can use cleanly:

prepped =
(baseDiscreteLHS[[1]]) //. {
e[1] -> e1,
e[2] -> e2,
delta["y"] -> dy,
delta["x"] -> dx,
UnitStep[r^2 - x^2 - y^2] -> usxy,
UnitStep[r^2 - x^2 - (-dy + y)^2] -> usxdy,
UnitStep[r^2 - (dx + x)^2 - y^2] -> usdxy,
UnitStep[r^2 - (dx + x)^2 - (-dy + y)^2] -> usdxdy,
UnitStep[r^2 - (-dx + x)^2 - y^2] -> usdxmdy
} // FullSimplify;


Now I'll use this stuff to pull the FD weights just for e1 from the first of the two equations provided:

weightsE1 =
{
fdmWeightDataND[prepped[[1]], e1, {x, y}, {dx, dy}],
fdmWeightDataND[prepped[[2]], e1, {x, y}, {dx, dy}],
fdmWeightDataND[prepped[[3]], e1, {x, y}, {dx, dy}],
fdmWeightDataND[prepped[[4]], e1, {x, y}, {dx, dy}],
fdmWeightDataND[
Simplify[Collect[prepped[[6]] // Expand, dx][[1]]],
e1,
{x, y},
{dx, dy}
],
fdmWeightDataND[
Simplify[Collect[prepped[[8]] // Expand, dx][[1]]],
e1,
{x, y},
{dx, dy}
]
};


Then I'll build all the matrices, add them together, and numericize everything:

mate1 =
Block[
{
r = .8, glass = (1.45)^2, air = 1., k = (2 \[Pi])/1.55,
sp1, spRules, points = 25,
dx, dy
},
dx = 1./points; dy = 1/points;
sp1 = FDM2D[#, {dx, dy}, {points, points}] & /@ weightsE1 // Total;
spRules =
# ->
With[{x = #[[1]]*dx, y = #[[2]]*dy},
#2 /.
{
usxy -> UnitStep[r^2 - x^2 - y^2],
usxdy -> UnitStep[r^2 - x^2 - (-dy + y)^2],
usdxy -> UnitStep[r^2 - (dx + x)^2 - y^2],
usdxdy -> UnitStep[r^2 - (dx + x)^2 - (-dy + y)^2],
usdxmdy -> UnitStep[r^2 - (-dx + x)^2 - y^2]
}
] &,
{
sp1["NonzeroPositions"],
sp1["NonzeroValues"]
}
];
SparseArray[spRules, Dimensions[sp1]]
]


Finally we can solve this using your same Eigensystem code:

Block[{r = .8, glass = (1.45)^2, air = 1., k = (2 \[Pi])/1.55, points = 25},
{val, vec} =
Eigensystem[mate1, -6, Method -> {"Arnoldi", "Shift" -> k^2 glass}
];
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All,
{mat, val}
]
] // AbsoluteTiming


The results are meaningless since I just randomly took part of the problem, but this entire procedure was fast and not memory intensive at all.

• Well, I should say the method is a bit too limited if it can only handle 1D case. Actually the memory issue mainly comes up in high dimensional cases. Commented Jun 19, 2019 at 10:27
• @xzczd my sense is the 1D case can be generalized out if you use a similar preprocessing approach to get the nD forms or if you can reduce the nD case to a tensor product of 1D cases with it. That’ll take work I didn’t have time for though. Commented Jun 19, 2019 at 15:46
• @xzczd I partially handled the nD case Commented Jun 27, 2019 at 19:19
• This code when points = 100 works 10 times longer than the code by @xzczd. And there is no result. What is the benefit here? Commented Jun 29, 2019 at 6:51
• (+51) 千金市骨 :) Commented Jun 29, 2019 at 9:54

To compare code @xzczd and @b3m2a1, I took a simplified version of the 2D operator from the topic here and 300x300 cells. I simplified code @xzczd a little, so it works 3 times faster. Code @b3m2a1, on the contrary, I complicated a little, so that it works 30 times faster. As a result, they now work equally fast and produce the same output.

Code @xzczd

r = .8; glass = (1.45)^2; air = 1.; k = (2 \[Pi])/1.55;
n2[x_, y_] := If[x^2 + y^2 <= r^2, glass, air];
ClearAll[fw, bw, delta]
SetAttributes[#, HoldAll] & /@ {fw, bw};

fw@D[expr_, x_] :=
With[{\[CapitalDelta] = delta@ToString@x},
Subtract @@ (expr /. {{x -> x + \[CapitalDelta]}, {x ->
x}})/\[CapitalDelta]]
bw@D[expr_, x_] :=
With[{\[CapitalDelta] = delta@ToString@x},
Subtract @@ (expr /. {{x -> x}, {x ->
x - \[CapitalDelta]}})/\[CapitalDelta]]

delta[a_ + b_] := delta@a + delta@b
delta[k_. delta[_]] := 0
grad = Function[expr, fw@D[expr, #] & /@ {x, y}];
div = Function[expr, Total@MapThread[bw@D@## &, {expr, {x, y}}]];
elst = e[#][x, y] & /@ Range[1];
discretelhs =
With[{n2 = n2[x, y]}, div@grad@elst + (k^2*n2 ) elst] //
FullSimplify;
points = 302; L = 2; domain = {-L, L}; grid =
Array[# &, points, domain];
del = #[[2 ;; -2]] &;
delta["x"] = delta["y"] = -Subtract @@ domain/(points - 1);
vare = Outer[e[#][#2, #3] &, Range@1, del@grid, del@grid, 1] //
Flatten;
discretelhslst =
Flatten@Transpose[
Module[{e}, e[i_][L | -L, _] = e[i_][_, L | -L] = 0.;
Table[discretelhs, {x, del@grid}, {y, del@grid}]], {2, 3,
1}]; // AbsoluteTiming
{barray2, marray2} =
CoefficientArrays[discretelhslst, vare]; // AbsoluteTiming
{val, vec} =
Eigensystem[marray2, -4,
Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All,
ColorFunction -> "AvocadoColors"] &, {mat, val}]


Code @b3m2a1

fdmWeightDataND~SetAttributes~HoldAll;
fdmWeightDataND[a_, dep_Symbol, indeps : {__Symbol},
hs : {__Symbol}] :=
Block[{dep},
Block[indeps,
Block[hs,
Module[{weights, order, denomWeight, denom, num, powed, pows,
minpow, pts, x}, denom = Denominator[a];
(*assumes the mesh spacing is in there at different powers*)
order = Exponent[denom, #] & /@ hs;
(*pull the power on the mesh spacing*)
denomWeight = If[order > 0, Coefficient[denom, hs^order], denom];
num = Numerator[a];
{powed, pows} =
Reap[num /.
dep[args___]?(Not@*FreeQ[Alternatives @@ indeps]) :>
dep @@ Sow[MapThread[Coefficient[#, #2] &, {{args}, hs}]]];
minpow = Min[pows];
weights =
Map[Coefficient[powed /. (dep @@ #) -> x, x] &, pows[[1]]]];
<|"Weights" -> weights, "Orders" -> order|>]]]];

FDM2D[weightData_, {hx_, hy_}, {nptsx_, nptsy_}] :=
Module[{xps, yps, weights = weightData["Weights"],
orders = weightData["Orders"]},
xps = Length@DeleteDuplicates@Keys[weights][[All, 1]];
yps = Length@Sort@DeleteDuplicates@Keys[weights][[All, 2]];
SparseArray[
Flatten@KeyValueMap[
With[{ix = #[[1]], iy = #[[2]], v = #2},
Map[Band[
If[ix >= 0, {1, 1 + ix}, {1 - ix, 1}] + {#, #} +
If[iy >= 0, {0, iy*nptsx}, {-iy*nptsx, 0}], {nptsx*nptsy,
nptsx*nptsx}, {nptsx, nptsx}] -> v &,
Range[nptsx - Abs[ix]] - 1]] &, weights], {nptsx*nptsy,
nptsx*nptsx}]/(Times @@ ({hx, hy}^orders))]

ClearAll[fw, bw, delta]
SetAttributes[#, HoldAll] & /@ {fw, bw};

fw@D[expr_, x_] :=
With[{\[CapitalDelta] = delta@ToString@x},
Subtract @@ (expr /. {{x -> x + \[CapitalDelta]}, {x ->
x}})/\[CapitalDelta]]
bw@D[expr_, x_] :=
With[{\[CapitalDelta] = delta@ToString@x},
Subtract @@ (expr /. {{x -> x}, {x ->
x - \[CapitalDelta]}})/\[CapitalDelta]]

delta[a_ + b_] := delta@a + delta@b
delta[k_. delta[_]] := 0
grad = Function[expr, fw@D[expr, #] & /@ {x, y}];
div = Function[expr, Total@MapThread[bw@D@## &, {expr, {x, y}}]];
elst = e[#][x, y] & /@ Range[1];
prepped = (baseDiscreteLHS[[1]]) //. {e[1] -> e1, delta["y"] -> dy,
delta["x"] -> dx};
weightsE1 = {fdmWeightDataND[prepped, e1, {x, y}, {dx, dy}]};
mate1 = Block[{sp1, spRules, points = 300, dx, dy},
dx = 4./(points + 1); dy = 4./(points + 1);
sp1 = FDM2D[#, {dx, dy}, {points, points}] & /@ weightsE1 // Total;
spRules =
MapThread[# -> With[{x = #[[1]]*dx, y = #[[2]]*dy}, #2] &, {sp1[
"NonzeroPositions"], sp1["NonzeroValues"]}];
SparseArray[spRules, Dimensions[sp1]]]; // AbsoluteTiming
r = .8; glass = (1.45)^2; air = 1.; k = (2 \[Pi])/1.55;
n2[x_, y_] := If[x^2 + y^2 <= r^2, glass, air]; discretelhs =
With[{n2 = n2[x, y]}, (k^2*n2) elst];
points1 = 302; L = 2; domain = {-L, L}; grid =
Array[# &, points1, domain];
del = #[[2 ;; -2]] &;
delta["x"] = delta["y"] = -Subtract @@ domain/(points1 - 1);
vare = Outer[e[#][#2, #3] &, Range@1, del@grid, del@grid, 1] //
Flatten;
discretelhslst =
Flatten@Transpose[
Module[{e}, e[i_][L | -L, _] = e[i_][_, L | -L] = 0.;
Table[discretelhs, {x, del@grid}, {y, del@grid}]], {2, 3, 1}];
{barray2, marray2} =
CoefficientArrays[discretelhslst, vare]; // AbsoluteTiming
points = 300; {val, vec} =
Eigensystem[mate1 + marray2, -4,
Method -> {"Arnoldi", "Shift" -> k^2*glass}]; // AbsoluteTiming
mat = ArrayReshape[#, {2, points, points}] & /@ vec;
ArrayPlot[#[[1]], PlotLabel -> Sqrt[#2], PlotRange -> All,
ColorFunction -> "AvocadoColors"] &, {mat, val}]


Finally, compare with my method that I suggested for this problem here using FEM:

<< NDSolveFEM
r = 0.8; ne = 4;
reg = ImplicitRegion[-2 <= x <= 2 && -2 <= y <= 2, {x, y}];
mesh = ToElementMesh[reg,
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.0004 (1 + 9 Norm[Mean[vertices]])]];

glass = 1.45; air = 1.; k0 = (2 \[Pi])/1.55; b = I*Sqrt[glass]*k0*1.37;
n[R_] := If[R <= r, glass^2, air^2]*k0^2

helm = -Laplacian[u[x, y], {x, y}] - (b^2 + n[Sqrt[x^2 + y^2]])*
u[x, y];
boundary = DirichletCondition[u[x, y] == 0, True];

{vals, funs} =
NDEigensystem[{helm, boundary}, u[x, y], {x, y} \[Element] mesh,
ne]; // AbsoluteTiming
Table[DensityPlot[Im[funs[[i]]], {x, y} \[Element] mesh,
PlotRange -> All,
PlotLabel -> Im[Sqrt[vals[[i]] + b^2]], Mesh -> None,
ColorFunction -> "AvocadoColors", Frame -> True, Axes -> False,
FrameTicks -> None], {i, Length[vals]}]


Comparing the data in Figure 1-3, we find that there is an eigenvalue match in methods FDM and FEM. The small difference can be explained by the difference in the number of cells.

• Very nice! What did you have to patch in my code, if you don't mind my asking? I should integrate the changes and quash whatever bugs you found. Commented Jul 1, 2019 at 7:35
• @b3m2a1 I used your code to calculate the laplacian matrix. I used a different code to calculate the diagonal matrix, since your code hangs on a diagonal matrix. Commented Jul 1, 2019 at 8:16

I'd like to add my incomplete answer. I call it incomplete because a general-purposed function hasn't been created (yet), but I do find a general strategy for handling the problem.

The idea is simple: generate coefficient array of general difference formula with CoefficientArray first, then generate the complete array using the general coefficients.

I'll solve the 2D problem in the question for illustration. delta, discretelhs etc. are the same as in the question.

sfx = Rescale[#, {x, x + delta@"x"}, {nx, nx + 1}] &;
sfy = Rescale[#, {y, y + delta@"y"}, {ny, ny + 1}] &;
rule = e[i_][x_, y_] :> {i, sfx@x // Simplify, sfy@y // Simplify};
varegeneral = Union@Cases[discretelhs, e[_][__], Infinity];
index = {i, nx, ny, ##} & @@@ (varegeneral /. rule);
coelst = Last@CoefficientArrays[discretelhs, varegeneral];
marrayfast =
Flatten[SparseArray[
Flatten[#,
3] &@(Unevaluated@
Compile[{},
Table[term, {i, 2}, {nx, 2, pts - 1}, {ny, 2, pts - 1}]] /. {term -> index,
pts -> points})[] -> (Flatten[#,
3] &@(Unevaluated@
Compile[{}, Table[lst[[i]], {i, 2}, {x, dgrid}, {y, dgrid}]] /. {dgrid ->
del@grid,
lst -> (Normal@coelst /. {HoldPattern@glass -> glass,
HoldPattern@air -> air})})[]), {2, points, points, 2, points,
points}][[All, 2 ;; -2, 2 ;; -2, All, 2 ;; -2, 2 ;; -2]], {{1, 2, 3}, {4, 5,
6}}]; // AbsoluteTiming

MaxMemoryUsed[]/1024^2. MB

marrayfast == marray2
(* True *)

{val, vec} =
Eigensystem[marrayfast, -6,
Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming

MaxMemoryUsed[]/1024^2. MB

mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
MapThread[ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All,
ColorFunction -> "AvocadoColors"] &, {mat, val}]


When points = 500, it only takes 5.63068 seconds to generate marrayfast, while 53.4027 + 19.1697 to generate marray2.

Memory usage is also reduced, but the real bottleneck seems to be Eigensystem, so perhaps it's not a big deal.

• My guess is you can generalize out your marrayfast using some type of KroneckerProduct. Unless that is where you are also choosing to deal with that pesky UnitStep` term. Commented Jun 29, 2019 at 16:41