Defining a big matrix for iterative algorithm

Question

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?

Answer 1

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...

Answer 2

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.

Answer 3

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