Defining a big matrix for iterative algorithm

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

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.

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