# writing linear system in matrix form

I have the following set of linear equations

$$a_{m+1,n}+a_{m-1,n}+m(a_{m,n+1}+a_{m,n-1})+(m^{2}+n^{2})a_{m,n}=f_{m,n}$$

Here $m$ and $n$ run from 1 to $N$, so there are $N^2$ equations for the unknowns $a_{m,n}$. Fixed boundary conditions are used such that $a_{N+1,n}=0$ etc.

I would like to solve this system numerically, and for that I need to first write it in matrix form $M a = f$. How can I programatically generate the coefficient matrix M?

EDIT Trying to follow @Wolfgang's suggestion

Ok, so lets take the simplest case of N=2. Writing the equations by hand and extracting the coefficients results in the matrix equation

$$\left[\begin{array}{cccc} 2 & 1 & 1 & 0\\ 1 & 5 & 0 & 1\\ 1 & 0 & 5 & 2\\ 0 & 1 & 2 & 8 \end{array}\right]\left[\begin{array}{c} a_{1,1}\\ a_{1,2}\\ a_{2,1}\\ a_{2,2} \end{array}\right]=\left[\begin{array}{c} f_{1,1}\\ f_{1,2}\\ f_{2,1}\\ f_{2,2} \end{array}\right]$$

Now if I follow @Wolfgang's suggestion and replace $(m,n)$ by $i$ I get a different matrix

$$\left[\begin{array}{cccc} 5 & 1 & 2 & 0\\ 1 & 5 & 1 & 1\\ 2 & 1 & 8 & 1\\ 0 & 1 & 1 & 10 \end{array}\right]\left[\begin{array}{c} a_{1}\\ a_{2}\\ a_{3}\\ a_{4} \end{array}\right]=\left[\begin{array}{c} f_{1}\\ f_{2}\\ f_{3}\\ f_{4} \end{array}\right]$$

• Minor nitpick: I would like to solve this system numerically, and for that I need to first write it in matrix form -- no, you don't need to do it. There are very popular matrix-free methods, which are often the best choice for large-scale systems. – Federico Poloni Oct 5 '14 at 13:45
• I've edited my answer in response to Andy Terrel's comment as pointed out below my answer. – Wolfgang Bangerth Oct 6 '14 at 15:27

You need to re-enumerate your degrees of freedom. Right now, you have two indices $1\le m,n\le N$. These span a rectangular array. Enumerate them from $i=0$ to $N^2-1$ (it turns out that zero-based indexing is so much simpler), for example by defining $i(m,n)=(m-1)+N*(n-1)$. In your case, this would corresponding to row-wise numbering. You can recover the original indices by observing that $m=(i\%N)+1$ and $n=(i/N)+1$ where the two operations are done in integer arithmetic. Then you can substitute these definitions in your formula everywhere.
Specifically, if you write your linear system in terms of $a_i$, not $a_{m,n}$, then you'll get this set of equations: $$a_{i+1}+a_{i-1}+(i\%N+1)(a_{i+N}+a_{i-N})+((i\%N+1)^2+(i/N+1)^2)a_i = f_i$$ where I've again assumed that you do the operations $i\%N$ and $i/N$ in integer arithmetic.
This is the equation you have for $i=0,\ldots,N^2$. Now write each of these $N^2$ equations down one by one and you get the corresponding matrix entries.
• @Andrei, Wolfgang has written the general form of the $i^{\rm th}$ row of the matrix. So if you take $i$ from 1 to $N^2$, you get the coefficients of the matrix on each row. – Bill Barth Oct 3 '14 at 14:08
• Not quite. The entries in $M_{1,2},M_{2,1},M_{3,4},M_{4,3}$ should be 1, if I read your formula correctly. – Wolfgang Bangerth Oct 3 '14 at 14:29
• If you define $i$ as $m+Nn$ then it runs from $N+1$ to $N(N+1)$ as $m$ and $n$ run from 1 to $N$, and not $1$ to $N^2$. Furthermore, I still cannot recover the matrix $M$ obtained by manual substitution. – Andrei Oct 4 '14 at 11:18