Implementing routine for $-\nabla\cdot (k(x,y) \nabla u)=f$ in Matlab

Question

I am solving the Poisson Equation for 2D given by the following expression: $$-\nabla\cdot (k(x,y) \nabla u)=f$$ in a rectangle with Dirichlet conditions on the boundary using Matlab. In principle I have made a code to solve the same problem, but in 1D and the code is as follows:

function [U,x]=p1d(n,alpha,beta,a,f)%alpha,beta, boundary conditions, here a is the k
   x=linspace(0,1,n);x=x(:);
   h=x(2)-x(1);
   pm=(x(1:end-1)+x(2:end))/2;
   K=diag(a(pm(1:end-1))+a(pm(2:end)))-diag(a(pm(2:end-1)),-1)-diag(a(pm(2:end-1)),1);
  F=h.^2*f(x(2:n-1));
  F(1)=F(1)+a(pm(1))*alpha;
  F(n-2)=F(n-2)+a(pm(end))*beta;
  U=K\F;
  U=[alpha;U;beta];
 endfunction

The problem I have to transfer the same ideas to 2D is the construction of the $K$ matrix. I understand how to add the conditions to the boundary, but the problem is that in the 2D case the position of each node is given by a strange shape, something like $l=(i-1)+(j-2)*(n-2)$.

By discretization, which is the following, using finite differences:

$-\Big(\dfrac{k_{i,j+\frac{1}{2}} u_{i,j+1}-\big(k_{i,j+\frac{1}{2}}+k_{i,j-\frac{1}{2}}\big)u_{i,j}+k_{i,j-\frac{1}{2}}u_{i,j-1}}{h^2}\Big)- \Big(\dfrac{k_{i+\frac{1}{2},j} u_{i+1,j}-\big(k_{i+\frac{1}{2},j}+k_{i-\frac{1}{2},j}\big)u_{i,j}+k_{i-\frac{1}{2},j}u_{i-1,j}}{h^2}\Big)=f_{i,j}$

you can see that the matrix is going to be diagonal by blocks, so the construction causes me problem.

Well what I started doing is the basics:

function U=poisson2d(x,y,k,f,g)% g condition to the border
n=size(x,1); % Mesh size
h=x(1,2)-x(1,1); %step size
K=sparse((n-2)*(n-2),(n-2)*(n-2)); %(this is the matrix causing me problem to assemble)
B=sparse((n-2)*(n-2),1);U=zeros(n); %initialize right side and solution
pmx=(x(:,1:end-1)+x(:,2:end))/2;%this is for the midpoints
pmy=(y(1:end-1,:)+x(2:end,:))/2;

Adding the input from the source function and the conditions to the boundary is no problem for me. But I can't implement the part of building K, I have a lot of trouble with the construction because of the position issue.

I hope you can help me, thank you in advance.

Edit

According to the discretization for the case of a 3x3 mesh, we would have the following coefficient matrix (not counting the boundary conditions). \begin{equation} \begin{pmatrix} \gamma & -\Gamma & 0 & -\Omega & 0 &0& 0& 0& 0\\ -\alpha & \gamma & -\Gamma& 0 &-\Omega& 0 &0& 0 &0\\ 0 & -\alpha & \gamma &0&0&-\Omega&0&0&0\\ -\beta &0 &0&\gamma&-\Gamma&0&-\Omega &0&0\\ 0&-\beta&0&-\alpha&\gamma&-\Gamma&0&-\Omega&0\\ 0&0&-\beta&0&-\alpha&\gamma&0&0&-\Omega\\ 0&0&0&-\beta&0&0&\gamma&\Gamma&0\\ 0&0&0&0&-\beta&0&-\alpha&\gamma&\Gamma\\ 0&0&0&0&0&-\beta&0&-\alpha&\gamma \end{pmatrix} \end{equation}

Where: $\gamma=\big(k_{i,j+\frac{1}{2}}+k_{i,j-\frac{1}{2}}+k_{i+\frac{1}{2},j}+k_{i-\frac{1}{2},j}\big)$

$\alpha=k_{i,j-\frac{1}{2}}$, $\Gamma =k_{i,j+\frac{1}{2}}$, $\beta=k_{i-\frac{1}{2},j}$, $\Omega=k_{i+\frac{1}{2},j}$

But why must the matrix be symmetric?

Answer 1

This is a (FDM) supplement to VoB's answer. You could write an extremely vectorized and optimized solver for your problem, but that is not a good first idea. Writing a for loop is easier, and once you identify the opportunities for optimization, you can implement them. Here is how I would go about it (in pseudo-code):

Assume $0\leq i\leq n$ and $0\leq j\leq m$, and let $x_{i,j}$ be the grid points

For each grid point x(i,j)
    Map the coordinate (i,j) to an index l: l = i + (j-1)*n (call this function coord2ind) 
    Populate the l-th row of the matrix K:
        If x(i,j) is on a Dirichlet boundary
            K(l,l) = 1, the rest of the row is zero
            Move to the next grid point
        If x(i,j) is on a Neumann boundary
            Set the row appropriately as the discretization dictates
            Move to the next grid point
        Otherwise
            Map the coordinate (i,j+1) to the index m using coord2ind
            Set K(l,m) = k(i,j+1/2)/h^2  --- needs to be computed from the function k
            ---- Similarly set the matrix entries

The right handside vector can be simultaneously populated. Let me know if you have further questions.

Answer 2

Probably not what OP was waiting for, but I think it could be pretty instructive and useful.

FEM codes use a much different approach to build the so-called stiffness matrix. In practice, they loop over elements and compute for each element small matrices (in your case if you use linear elements 3by3) which are distributed to the right entries of the global stiffness matrix. Such a process is called assembly, and can be applied to general geometries, not only on rectangles like you have to do. Moreover, you avoid all that nasty indexing problems.

If you just need to solve your problem, you can have a look at step-5 from deal.II tutorials, a C++ finite element library which has awesome tutorials with lots of details.

In particular, there you have exactly your problem, modulo $f$ and the coefficient $k(x)$.