Memory Issue with Large Matrix Iterations

Question

I am working on a Jacobi/SOR problem with matrix iterations. I have a matrix of size n = 40, 80, however, I cannot get my iterations past 14 before my program slows down dramatically. Is there a way to speed this up?
n = 40;

(*Matrix information*)
Subscript[A, 2] = 
  Total[{DiagonalMatrix[Array[-1 &, n - 1], -1], 
    DiagonalMatrix[Array[2 &, n]], 
    DiagonalMatrix[Array[-1 &, n - 1], 1]}];
b = ConstantArray[0, n];
k = 1;
h = (1/(n + 1));
While[k < n + 1, b[[k]] = h^2 Sin[3 [Pi] (h*k)]; k++]

L = Total[{DiagonalMatrix[Array[-1 &, n - 1], -1], 
    DiagonalMatrix[Array[0 &, n]], 
    DiagonalMatrix[Array[0 &, n - 1], 1]}];
U  =  Total[{DiagonalMatrix[Array[0 &, n - 1], -1], 
    DiagonalMatrix[Array[0 &, n]], 
    DiagonalMatrix[Array[-1 &, n - 1], 1]}];
d  =  Total[{DiagonalMatrix[Array[0 &, n - 1], -1], 
    DiagonalMatrix[Array[2 &, n]], 
    DiagonalMatrix[Array[0 &, n - 1], 1]}];
Subscript[C, J] = Dot[-Inverse[d], L + U];
No = Norm[Subscript[C, J]];

Timing[
 
 (*The Jacobi iteration loop*)
 x[0] = ConstantArray[0, n];
 x[n_] := x[n] = Dot[Subscript[C, J], x[n - 1]] + Dot[Inverse[d], b];
k = 1;

While[k < 6000, 
 Print["At k = ", k, " the approx is ", N[x[k]], ". Error is: ", 
  N[(No/(1 - No))*Norm[x[k] - x[k - 1]]]];
 If[(No/(1 - No))*Norm[x[k] - x[k - 1]] < (10)^-3, 
  Print["Success @ ", k]; Break[], k++]
 ]
 ]
```

Cantorexample · Accepted Answer

Here is what was wrong and how to fix the issue:
Problem:
Mathematica was storing values symbolically, which means it was storing functions to be evaluated that were of lengths 3000+ in a single array cell. This is unlike other programs like Python and Matlab, which do things numerically automatically. The symbolic evaluation in Mathematica is why the matrix iterations slowed down at a certain point; the concatenation of freakishly long functions to be evaluated.
Solution:
The solution to this is exceptionally simple, and, admittedly, I should have caught it much sooner. Instead of:
 (*The Jacobi iteration loop*)
 x[0] = ConstantArray[0, n];
 x[n_] := x[n] = Dot[Subscript[C, J], x[n - 1]] + Dot[Inverse[d], b];

It should be:
x[0] = ConstantArray[0, n];
 x[n_] := x[n] = N[Dot[Subscript[C, J], x[n - 1]]] + N[Dot[Inverse[d], b]];

Simply adding the N converts to numerical storing. With this, Mathematica is back to high speed.
Conclusion:
ALWAYS STORE YOUR VALUES NUMERICALLY WHEN DOING NUMERICAL METHODS!!!

Memory Issue with Large Matrix Iterations

One Answer

Add your own answers!

Ask a Question