# How to construct a time-dependent matrix quickly?

In the process of discretization of a 4D PDE, I need to construct a final sparse matrix $$B$$, which is very large (I denoted by $$\text{size}$$ here) and time-dependent, viz, some of its entries change with $$t$$. To move into the 4D discretization matrices, I use the command KroneckerProduct[].

A sample case is as follows:

ClearAll["Global*"];
SeedRandom;

m1 = 8; m2 = 8; m3 = 8; m4 = 8; size = m1*m2*m3*m4;
gamma = 0.3; etad = 0.7/100; etaf = 1.2/100;
rosv = -0.4; rosd = -0.15; rovd = 0.3; rosf = -0.15; rovf = 0.3; rodf = 0.25;
e = 100.; κ = 0.5; η = 0.1; lambdad = 0.01; lambdaf = 0.05;

c1 = 0.074; c2 = 0.014; c3 = 2.10; thetad[t] := c1 - c2*Exp[-c3*t]
f1 = 1.0; f2 = 0.5; f3 = 0.5; thetaf[t] := f1 - f2*Exp[-f3*t]

nx = Sort@RandomReal[{0, 100}, m1];
ny = Sort@RandomReal[{0, 10}, m2];
nz = Sort@RandomReal[{-1, 1}, m3];
nw = Sort@RandomReal[{-1, 1}, m4];
origrid = Flatten[Outer[List, nx, ny, nz, nw], 3];

Idx = SparseArray[{{i_, i_} -> 1.}, {m1, m1}, 0]; Idy =
SparseArray[{{i_, i_} -> 1.}, {m2, m2}, 0];
Idz = SparseArray[{{i_, i_} -> 1.}, {m3, m3}, 0]; Idw =
SparseArray[{{i_, i_} -> 1.}, {m4, m4}, 0];
DS = KroneckerProduct[(SparseArray@DiagonalMatrix@nx), Idy, Idz, Idw];
DV = KroneckerProduct[Idx, (SparseArray@DiagonalMatrix@ny), Idz, Idw];
DRD = KroneckerProduct[Idx, Idy, (SparseArray@DiagonalMatrix@nz), Idw];
DRF = KroneckerProduct[Idx, Idy, Idz, (SparseArray@DiagonalMatrix@nw)];
Id = KroneckerProduct[Idx, Idy, Idz, Idw];

dudx = SparseArray[{Band[{1, 1}] -> 1., Band[{2, 1}] -> 1.5,
Band[{1, 2}] -> -0.3}, {m1, m1}];
d2udx2 = SparseArray[{Band[{1, 1}] -> 2., Band[{2, 1}] -> 1.,
Band[{1, 2}] -> 0.6, Band[{3, 1}] -> 1.}, {m1, m1}];
dudy = SparseArray[{Band[{1, 1}] -> 0.75, Band[{2, 1}] -> 5.,
Band[{1, 2}] -> 2.}, {m2, m2}];
d2udy2 = SparseArray[{Band[{1, 1}] -> 3.5, Band[{2, 1}] -> 1.,
Band[{1, 2}] -> 2., Band[{3, 1}] -> 0.1}, {m2, m2}];
dudz = SparseArray[{Band[{1, 1}] -> 1.2, Band[{2, 1}] -> 1.,
Band[{1, 2}] -> 2.}, {m3, m3}];
d2udz2 = SparseArray[{Band[{1, 1}] -> 1.5, Band[{2, 1}] -> 1.,
Band[{1, 2}] -> 1., Band[{3, 1}] -> 0.3}, {m3, m3}];
dudw = SparseArray[{Band[{1, 1}] -> 5., Band[{2, 1}] -> 2.,
Band[{1, 2}] -> 1.}, {m4, m4}];
d2udw2 = SparseArray[{Band[{1, 1}] -> 3., Band[{2, 1}] -> 0.3,
Band[{1, 2}] -> 2., Band[{3, 1}] -> 0.2}, {m4, m4}];

B0 = SparseArray[
+(1/2 (DS^2).(DV)).KroneckerProduct[d2udx2, Idy, Idz, Idw]
+ (1/2 gamma^2 DV).KroneckerProduct[Idx, d2udy2, Idz, Idw]
+ (1/2 etad^2)*KroneckerProduct[Idx, Idy, d2udz2, Idw]
+ (1/2 etaf^2)*KroneckerProduct[Idx, Idy, Idz, d2udw2]
+ ((rosv*gamma)*(DS.DV)).KroneckerProduct[dudx, dudy, Idz, Idw]
Idw]
dudz, Idw]
+ ((rosf*etaf)*(DS.Sqrt[DV])).KroneckerProduct[dudx, Idy, Idz,
dudw]
+ ((rovf*gamma*etaf)*(Sqrt[DV])).KroneckerProduct[Idx, dudy, Idz,
dudw]
+ ((DRD - DRF).DS).KroneckerProduct[dudx, Idy, Idz, Idw]
+ (κ (η*Id - DV)).KroneckerProduct[Idx, dudy, Idz,
Idw]
- DRD.Id
]; // AbsoluteTiming
matcoef1 = (SparseArray[{{i_, i_} -> c1 - c2*Exp[-c3*t]}, {size,
size}, 0.]);
matcoef2 = (SparseArray[{{i_, i_} -> f1 - f2*Exp[-f3*t]}, {size,
size}, 0.]);
B[t_] := B0 + (lambdad (matcoef1 - DRD)).KroneckerProduct[Idx, Idy,
dudz, Idw] + ((lambdaf (matcoef2 - DRF)) -
rosf*etaf*Sqrt[DV]).KroneckerProduct[Idx, Idy, Idz,
dudw]; // AbsoluteTiming

B = B[t]; // AbsoluteTiming
MatrixPlot[B]


I know that it is correct, but is not at all efficient. For example, for the above case, only constructing $$B$$ takes around 27 seconds in my laptop, while I should choose much more inputs for $$m_1,m_2,m_3,m_4$$ in my real calculations.

I think the most time-consuming case is where at the end I involve matcoef1 and matcoef2 into calculations. Please note that I need to save $$B$$ since I will impose the boundary conditions later (based on a combination of its rows).

Does anyone know how we can speed up the process of constructing sparse matrices involving time-dependent entries for large sizes? I also tried to use Module or Compile, but I failed!

• It is absolutely not clear to me why it should not be sufficient to compute B[t] only for numeric t. Using symbolic matrices of this size just feels wrong. – Henrik Schumacher Nov 16 '18 at 21:18
• Thanks for your care. You're right, it seems strange, but I don't have any other options. For example, consider that I should add a multiplier of rows 1000-1010 into the rows 2050-2060! I mean I impose the boundary conditions of the PDE. And then, I will get a final B[t] and give it to NDSolve. Because of this, I need to compute this fast. However, even if I do not impose the boundaries, NDSolve` will need B[t] symbolically to do the time stepping. Am I right? If you know any tips, please write. – Fazlollah Nov 16 '18 at 21:37