1
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Federico Poloni Oct 5 '14 at 13:45
  • $\begingroup$ I've edited my answer in response to Andy Terrel's comment as pointed out below my answer. $\endgroup$ – Wolfgang Bangerth Oct 6 '14 at 15:27
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I still don't see how I can get the matrix form of the linear system from this.. $\endgroup$ – Andrei Oct 3 '14 at 13:55
  • $\begingroup$ @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. $\endgroup$ – Bill Barth Oct 3 '14 at 14:08
  • $\begingroup$ I moved my comment to the original question, due to lack of space here. $\endgroup$ – Andrei Oct 3 '14 at 14:21
  • $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Oct 3 '14 at 14:29
  • $\begingroup$ 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. $\endgroup$ – Andrei Oct 4 '14 at 11:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.