Dear Wolfarm community, I am very new with Mathematica. I've written a code which works pretty well, but the big problem is that it is so slow that I cannot use it. I usually use Matlab for this type of coding. The equivalent Matlab version is thousand times faster. Here in the following I report a dummy (very much simplified) version of code that contains the main operations I need to perform. Maybe I make some basic mistakes in programming with Mathematica that dramatically slow down the execution. When the parameters n and m increase, the code becomes unusable. Note also that I tried to replace all the functions defined in the code with compiled functions with command Compile. I expected to see big improvements but unfortunately it was not so. Is there a more efficient way to write this code? Thank you very much! Enzo
Clear["Global`*"]
(*Differential equation*)
Eq = D1w[1, 1] + D2w[1, 1, 2]
D1w[1, 1] + D2w[1, 1, 2]
(*Define BSpline parameters*)
p = 2;
q = 2;
n = 5;
m = 5;
r = n + p + 1;
s = m + q + 1;
rint = r + 1 - 2 (p + 1);
sint = s + 1 - 2 (p + 1);
Uint = Range[0, 1, 1/(rint + 1)];
Vint = Range[0, 1, 1/(sint + 1)];
U = N[Join[ConstantArray[0, p], Uint , ConstantArray[1, p]]];
V = N[Join[ConstantArray[0, q], Vint , ConstantArray[1, q]]];
Uc = Table[(Sum[U[[h + k]], {k, 1, p}])/p, {h, 1, n + 1}];
Vc = Table[(Sum[V[[h + k]], {k, 1, q}])/q, {h, 1, n + 1}];
(*Define BSpline approximations (Discretization) *)
D1w11[ic_, jc_, i_, j_] :=
D[BSplineBasis[{p, U}, i - 1, u], {u, 1}] BSplineBasis[{q, V}, j - 1, v] W1[
i, j] //. {u -> Uc[[ic]], v -> Vc[[jc]]};
D2w112[ic_, jc_, i_, j_] :=
D[BSplineBasis[{p, U}, i - 1, u], {u, 1}] D[
BSplineBasis[{q, V}, j - 1, v], {v, 1}] W1[i, j] //. {u -> Uc[[ic]],
v -> Vc[[jc]]};
(*Substitute the above approximate derivatives into the original differential equation*)
subw11[ic_, jc_, i_, j_] := {D1w[1, 1] :> D1w11[ic, jc, i, j]};
subw112[ic_, jc_, i_, j_] := {D2w[1, 1, 2] :> D2w112[ic, jc, i, j]};
sub[ic_, jc_, i_, j_] := Join[subw11[ic, jc, i, j], subw112[ic, jc, i, j]];
EqDisc[ic_, jc_, i_, j_] := Eq //. sub[ic, jc, i, j];
(*Now collects terms in w1*)
KEq[ic_, jc_, i_, j_] := Coefficient[EqDisc [ic, jc, i, j], {W1[i, j]}];
KEqW1[ic_, jc_, i_, j_] := KEq[ic, jc, i, j][[1]];
(*Transform double -index formulation into single - index formulation*)
iMat = Table[Range[1, n + 1, 1], m + 1];
ivec = Flatten[ArrayReshape[iMat, {1, (n + 1) (m + 1)}]];
jMat = Transpose[Table[Range[1, m + 1, 1], n + 1]];
jvec = Flatten[ArrayReshape[jMat, {1, (n + 1) (m + 1)}]];
KKEqW1[r_, s_] :=
KEqW1 @@ List[ivec [[r]], jvec [[r]], ivec [[s]], jvec [[s]]];
(*Finally, create a matrix*)
K = Table[KKEqW1[r, s], {r, 1, (n + 1) (m + 1)}, {s, 1, (n + 1) (m + 1)}];