Constructing an arbitrary dimensional matrix with a certain rule

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

$$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\\$$

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!

• Is this homework? Why are are using Mathematica for your problem? – gwr Jan 29 '16 at 16:20
• It's not a homework per se. For our research, I am solving a book's problem (which uses MATLAB) one by one. And we decided to do everything in Mathematica, so we're sticking to mathematica. Now and then, I can't solve a problem. We are even worse in other programming languages, so we are sticking to mathematica. – Student Jan 29 '16 at 16:24
• Although your matrix is not sparse, perhaps you could still construct it using SparseArray and rules, then converting the SparseArray object to a Normal matrix. – MarcoB Jan 29 '16 at 16:24
• Trust me, I am trying for a long time and I will keep trying. But if someone comes along, and does it in 5 minutes, then I have nothing to complain. – Student Jan 29 '16 at 16:26
• As you should know, Mathematica list indexes start at one, not zero. At least you may have converted your expressions that way. – Dr. belisarius Jan 29 '16 at 16:37

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 !

• Okay, this one I just checked and it works!! I will understand the code now. Thanks Andre!! Yes, the other one was giving a wrong answer. Yours is correct. – Student Jan 29 '16 at 17:52

(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!

• +1. I though that indices for matrices in general start with 1? (Just to start one of those religious discussions... :) – gwr Jan 29 '16 at 16:45
• @gwr Eh, some people may disagree there: Edsger W. Dijkstra's article "Why numbering should start at zero". I have also heard the argument that, for a base $b$ the first $b^N$ non-negative integers are represented by exactly $N$ digits (including leading zeros) only if the numbering starts at $0$, so that's the convention that makes the most sense on computing systems. I, for one, don't care either way :-) – MarcoB Jan 29 '16 at 16:57
• That's the zeroest time I hear about this. ;-) – gwr Jan 29 '16 at 17:04
• @MarcoB: Actually there are probably some mistakes. First of all, I think in the definition of c_i that you wrote, it will be i=1 as the first condition instead of i=0. Also, for the [i,j] component, both are c[i]. Even though I corrected them, the problems are still there. For example, for n=1, it's a 2X2 matrix with M_{12}=-1/2 and M_{21}=1/2. Using the code you wrote, I am getting -1 and +1 respectively. – Student Jan 29 '16 at 17:50
• @Student Very possible. I suggest that you amend the definitions in your original question by taking into consideration that Mathematica's array indices start at $1$. That will make it much easier to troubleshoot my code. – MarcoB Jan 29 '16 at 18:02

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}]]}]


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

• I upvoted both the answers and accepted them. I will run, check and understand them now. But I am really grateful to see so many people gave some thoughts to my question. Thank you. – Student Jan 29 '16 at 17:24
• @Student. I fixed my solution. (It had the same problem as MarcoB's.) – march Jan 29 '16 at 18:32