Constructing an arbitrary dimensional matrix with a certain rule

Question

I want to construct the following matrix $M$ of dimension $(N+1)\times(N+1)$ with the following rule (for arbitrary N). Below, in the subscript, as usual, the 1st index gives the row and the second one the column. The subscript indices go from 0 to $N$.

\begin{equation} M_{00}=\frac{2N^2+1}{6},\, M_{NN}=-\frac{2N^2+1}{6}\\ M_{jj}=-\frac{x_j}{2(1-x_j^2)}\quad\mbox{for}\quad j=1\dots N-1\\ M_{ij}=\frac{c_i}{c_j}\frac{(-1)^{i+j}}{(x_i-x_j)}\quad\mbox{for}\quad i\neq j,\, i,j=0\dots N\\ \end{equation}

Above, $c_i=2$ for $i=0$ or $N$, and 1 otherwise. The $x_j$'s ($j=0\dots N$) are defined as $x_j=\cos(j\pi/N)$.

The above rule should unambiguously tell us all the matrix entries, but how to implement it in Mathematica is beyond me. A nice, understandable program will be very helpful. Thanks a lot in advance!

Answer 1

(Update: fixed incorrectly translated rules)

Remember that array indices start at 1 in Mathematica, so you have to adapt your definitions accordingly. Perhaps you could have started there...

You can quite literally translate your requirements into conditions and feed them to SparseArray, then use Normal to get a standard representation:

Clear[matrixgenerator]

matrixgenerator[n_Integer] :=
 Module[{x, c},
  x[j_] := Cos[(j - 1) Pi/n];
  c[i_] := If[i == 1 || i == n + 1, 2, 1];
  Normal@Simplify@SparseArray[{
      {1, 1} -> (2 n^2 + 1)/6,
      {n + 1, n + 1} -> -(2 n^2 + 1)/6,
      {j_, j_} /; 1 < j <= n -> -x[j]/2/(1 - x[j]^2),
      {i_, j_} /; i != j -> c[i] (-1)^(i + j)/c[j]/(x[i] - x[j])
      }]
 ]

matrixgenerator[5] // MatrixForm

SparseArray can be a very powerful tool even if you are not dealing with sparse arrays at all!

Answer 2

Here is a solution without the mess of reindexing:

NN = 5;

c[i_] = If[i == 0 || i == NN, 2, 1];
x[j_] = Cos[j Pi/NN];

result = Table[Switch[{i, j},
    {0, 0}, (2 NN^2 + 1)/6,
    {NN, NN}, -(2 NN^2 + 1)/6,
    {xy_, xy_}, -x[j] / (2 (1 - x[j]^2)),
    {_, _}, c[i] (-1)^(i + j)/(c[j] (x[i] - x[j]))],
   {i, 0, NN},
   {j, 0, NN}
   ];

result // Simplify // TableForm

Result different from the other answer ! To be verified !

Answer 3

Here is a quick Mathematica translation of Trefethen's cheb.m MATLAB routine for generating a Chebyshev spectral differentiation matrix:

chebmat[n_Integer?Positive, prec_: MachinePrecision] := 
 Module[{fac, xm},
        xm = ConstantArray[N[-Cos[π Range[0, n]/n], prec], n + 1];
        fac = Flatten[{2, PadRight[{}, n - 1, {-1, 1}], 2 - 4 Mod[n, 2]}];
        xm = Outer[Times, fac, 1/fac]/(xm - Transpose[xm] + IdentityMatrix[n + 1]);
        xm - DiagonalMatrix[Total[xm, {2}]]]

For instance, FullSimplify[chebmat[5, ∞]] quickly generates the matrix in andre's answer. However, a more appropriate demonstration of the use of the differentiation matrix would be to use it for what its name is implying. ;)

Here goes:

tst[x_] := Exp[x] Sin[5 x]
With[{n = 20},
     chebnodes = Cos[π Range[0, n]/n];
     diff = Reverse[chebmat[n].tst[chebnodes]]];

Plot[tst'[x], {x, -1, 1}, 
     Epilog -> {Directive[AbsolutePointSize[5], ColorData[97, 2]], 
                Point[Transpose[{-chebnodes, diff}]]}]

Answer 4

Alternatively to MarcoB's solution (but using their definitions), use a Table and a Which or a Piecewise. This is strictly for the purpose of showing more methods.

matrixgenerator[n_Integer] :=
 Module[{x, c}, x[i_] := Cos[i Pi/n];
   c[i_] := If[i == 0 || i == n, 2, 1];
   Table[
    Which[
     i == j == 0, (2 n^2 + 1)/6,
     i == j == n, -(2 n^2 + 1)/6,   
     i == j, -x[j]/(2 (1 - x[j]^2)),
     i != j, c[i] (-1)^(i + j)/c[j]/(x[i] - x[j])
     ],
    {i, 0, n}, {j, 0, n}]
   ] // Simplify

matrixgenerator[n_Integer] :=
 Module[{x, c}, x[i_] := Cos[i Pi/n];
   c[i_] := If[i == 0 || i == n, 2, 1];
   Table[
    Piecewise[{
      {(2 n^2 + 1)/6, i == j == 0},
      {-(2 n^2 + 1)/6, i == j == n},    
      {-x[j]/(2 (1 - x[j]^2)), i == j},
      {c[i] (-1)^(i + j)/c[j]/(x[i] - x[j]), i != j}
      }],
    {i, 0, n}, {j, 0, n}]
   ] // Simplify