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

1 Answer 1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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