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

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{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}$$$$

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?

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.

• Hello Abdullah, I have achieved the implementation. But I have a doubt a little bit more about that matrix K. This matrix must be symmetric?
– Haus
Jun 13 at 17:36
• Well, of course taking the fact that the k function must be symmetrical
– Haus
Jun 13 at 17:37
• It will not be completely symmetric. However, it should be symmetric after you eliminate the rows corresponding to the Dirichlet b.c.s. Jun 13 at 18:01
• If you pick $k(x,y) = 1$ (or any constant for that matter), the matrix you get seems to be correct. But depending on $k(x,y)$, the matrix may not be correct. Did you test your implementation with a manufactured solution to see if the numerical solution is close to the exact solution? Jun 13 at 20:10
• Compare these errors to errors you get when you pick $n=25$ and $n=100$. See if it is reducing optimally => halving the mesh size should (approximately) quarter the error. Jun 13 at 20:43

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

• Thank you very much, the main problem is my lack of knowledge in matlab, but as you mention, it comes very detailed. I will try to move it to Matlab.
– Haus
May 29 at 14:34
• In addition, what you mention is interesting about what works for domains that are not necessarily rectangular, I think it makes all the sense in the world, because staying only with rectangular domains would be very little, but it is my first course in relation to this, I hope at some point to deepen in that field of study
– Haus
May 29 at 14:37
• @Haus since this is for a course, I would suggest you to stick to the discretization you learned in class. If you use FEM but only covered FDM, you may miss important learning objectives. May 29 at 15:32
• @Haus If this is a homework, then of course you have to stick to what you learned, as Adbullah wrote. I assumed you were using finite elements since you used that question tag, but if you must use FDM, my answer is definitely not an answer at all.
– VoB
May 29 at 16:08
• Thanks to both. Yes, it is based on FDE and it is a task, only that I have been trying to do the assembly of the matrix for days and I can not do it
– Haus
May 29 at 17:20