# Defining a big matrix for iterative algorithm

I need to solve this system of equations with Seidel and Jacobi methods, so I tried to define a big matrix in Wolfram.

I have read documentation, and I wrote the following command that should help me, but it is now working:

Table[If[i == j, a], If[i == j + 1, 1], If[i + 1 == j, 1],
If[i + 2 == j, 1/b], {i, 100}, {j, 100}]


Can we work in wolfram with such big matrixes? If yes, how I should declare them?

And maybe someone can advise me where I can read about iterative algorithms on matrixes in Wolfram?

• If we are using decimal values rather than exact values, 100 x 100 is actually not at all large. – Daniel Lichtblau Oct 13 at 22:00

I'll illustrate with a smaller example, using the Jacobi method. Since the matrix is sparse it makes sense to use the SparseArray structure from the Wolfram Language.

a = b = 10.;
n = 16;
mat = SparseArray[{{i_, i_} -> a, {i_, j_} /; j == i - 1 :>
1., {i_, j_} /; j == i + 1 :> 1., {i_, j_} /; j == i + 2 :> 1/b}, n];
rhs = N@Range[n];


For the iterations we separate the (dominant) diagonal and the rest of the matrix.

diag = Normal[Diagonal[mat]];
rest = mat - DiagonalMatrix[diag];


We'll initialize a result and do a few Jacobi iterations.

x[0] = ConstantArray[0., n];
Do[x[j] = 1/diag*(rhs - rest.x[j - 1]), {j, 4}];


What is this approximated result now?

(* In[512]:= *)x[4]

(* Out[512]= {0.0809433, 0.1635632, 0.2460271, 0.328511, 0.4109949, \
0.4934788, 0.5759627, 0.6584466, 0.7409305, 0.8234144, 0.9059, \
0.988435, 1.071278, 1.15244, 1.2357, 1.47554} *)


Here is what a direct linear solver gives.

soln = LinearSolve[mat, rhs]

(* Out[518]= {0.0811406188159, 0.163937166886, 0.24656644955, \
0.329212627692, 0.411857099219, 0.494501743134, 0.577146369797, \
0.659790996402, 0.742435624897, 0.825080412827, 0.90772341803, \
0.990368288034, 1.07317118838, 1.15422513249, 1.23694695631, \
1.47630530437} *)


We can see they are already not too far apart.

Norm[x[4] - soln]

(* Out[522]= 0.00519376049093 *)


A better iteration would use e.g. FixedPoint and stop when the change is less than some epsilon.

By the way, this can be carried out in exact arithmetic as well. Expect large numerators and denominators...

• Thank you very much! – Yura Holubeu Oct 14 at 9:56

# New Method

n = 100;
method2[n_] := SparseArray[{
{i_, i_} -> a,
{i_, j_} /; (i == j + 1 || i == j - 1) -> 1,
{i_, j_} /; i == j - 2 -> 1/b},
{n, n}
]

method1[n].Array[x, n] // RepeatedTiming;
method2[n].Array[x, n] // RepeatedTiming;
{#[[;; , 1]], SameQ @@ #[[;; , 2]]} &@Out[{-2, -1}]


{{0.020, 0.0062}, True}

# Original Method

n = 10;
method1[n_] := Table[
Which[i == j, a, i == j + 1, 1, i + 1 == j, 1, i + 2 == j, 1/b, True, 0],
{i, n}, {j, n}
];
method1[n].Array[x, n] // TableForm

(* a x[1]+x[2]+x[3]/b
x[1]+a x[2]+x[3]+x[4]/b
x[2]+a x[3]+x[4]+x[5]/b
x[3]+a x[4]+x[5]+x[6]/b
x[4]+a x[5]+x[6]+x[7]/b
x[5]+a x[6]+x[7]+x[8]/b
x[6]+a x[7]+x[8]+x[9]/b
x[7]+a x[8]+x[9]+x[10]/b
x[8]+a x[9]+x[10]
x[9]+a x[10] *)


As for the resources on iterative methods, I'll keep looking.

• Thank you for your answer! But can this code add each element from each sum to a separate element of matrix? – Yura Holubeu Oct 13 at 20:30
• @YuraHolubeu Sorry, I don't understand the question. – That Gravity Guy Oct 13 at 21:21
• Thanks, I will look to this code as soon as possible :) – Yura Holubeu Oct 14 at 0:10

An alternative way to construct a SparseArray using Band:

sa[n_] := SparseArray[{Band[{1, 1}] -> a,
Band[{2, 1}] -> 1,
Band[{1, 2}] -> 1,
Band[{1, 3}] -> 1/b},
{n, n}]

Array[x, 10] # & /@ sa[10] // MatrixForm // TeXForm


$$\small\left( \begin{array}{cccccccccc} a x(1) & x(2) & \frac{x(3)}{b} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x(1) & a x(2) & x(3) & \frac{x(4)}{b} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & x(2) & a x(3) & x(4) & \frac{x(5)}{b} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & x(3) & a x(4) & x(5) & \frac{x(6)}{b} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & x(4) & a x(5) & x(6) & \frac{x(7)}{b} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & x(5) & a x(6) & x(7) & \frac{x(8)}{b} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & x(6) & a x(7) & x(8) & \frac{x(9)}{b} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & x(7) & a x(8) & x(9) & \frac{x(10)}{b} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & x(8) & a x(9) & x(10) \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & x(9) & a x(10) \\ \end{array} \right)$$

Note: Using Band is faster than method2 in ThatGravityGuy's answer:

n = 100;
sa[n].Array[x, n]; // RepeatedTiming// First


0.00051

method2[n].Array[x, n]; //RepeatedTiming// First


0.0071

• Thank you very much, this looks like an ideal answer! – Yura Holubeu Oct 14 at 9:56
• @YuraHolubeu, you are welcome. And welcome to mma.se. – kglr Oct 14 at 10:13