# Continuation of solution to $-\nabla\cdot (k(x,y)\nabla u)=f$

I'm trying to solve the following problem, I had previously opened another discussion for the implementation and well, it seems that it has turned out well, it can be found here.

I need to calculate the preconditioned conjugate gradient solution with incomplete Cholesky and Jacobi, for a mesh of n = 500,1000,5000 and 10,000 points. But these are a lot of points for a personal computer. As @Abdullah Ali Sivas comments:

You don't need the matrix, you can use the action of the matrix to solve the Poisson problem; especially if you implement stationary iterative methods like Jacobi or Gauss-Seidel methods.

I would like you to please help me to know how to obtain the solution to such a problem when the dimensions are gigantic, there are many points per side and my computer does not have such capacity.

Using Handle function:

function y=afun(x)
nx = 7;
N = nx*nx;
x1 = linspace(0,1,nx);
y1 = x1;
[X,Y] = ndgrid(x1,y1);
dx = x1(2)-x1(1);
k=@(x,y) 1+x.^2+y.^2;
isDirichlet = (X==0) | (X==1) | (Y==0) | (Y==1);
Ad(:,5) =  -k( X(:), Y(:)-dx/2 );
idx = find(isDirichlet(:));
n=sqrt(size(C,1));
l=B(1:end-1);s=U(2:end);
for i=n:n:(n-1)*n
l(i)=0;s(i)=0;
end
y=([0;l]+[C]+[s;0]+[diag(zeros(n));L(1:(n-1)*n)]+[UP(n+1:end);diag(zeros(n))]).*x;
end

• What is the form of k(x,y)? The problem is linear, and for some k(x,y) and some BC one could use FFT in one of the coordinates, and solve for one Fourier harmonic at a time; that would lead to a set of smaller linear problems (with complex coefficients). – Maxim Umansky Jun 21 at 4:37
• $k(x,y)$ is a function of a real and symmetric variable – Haus Jun 21 at 4:44
• What is symmetric variable? Can you write explicitly some possible k(x,y)? – Maxim Umansky Jun 21 at 5:41
• Of course, $k(x,y)=1+x^2+y^2$, symmetry as function. – Haus Jun 21 at 5:41
• Is the domain circular by any chance? Or rectangular? – Maxim Umansky Jun 21 at 14:10

According to your first question you have a sparse matrix. That means that you don't need to store all the $$n\times n$$ coefficients for your problem. A format that is useful for matrix-vector products is Compressed Sparse Row (CSR).

In the case of the Jacobi method you can rearrange

$$Ax = f\, ,$$

as

$$x^{k + 1} = D^{-1}(f - (L + U) x)\, ,$$

with $$D$$ the diagonal of $$A$$, and $$L + U = A - D$$.

The same treatment goes for Gauss-Seidel or Conjugate Gradient methods.

• But with this we would not need to calculate the coefficients of A? – Haus Jun 21 at 14:37
• @Haus Yes. And store it in a sparse matrix, or several, depending on the method. – nicoguaro Jun 21 at 14:40

First of all, I want to ask if you are sure that when you are saying $$n=10000$$ you mean $$10^4$$ points in one dimension? Because that would mean $$10^8$$ points over the domain, and that is too many points to solve such a simple problem. You would hit the minimum error before $$10^8$$ points, and you wouldn't gain anything for the extra work you do.

Secondly, you don't really need the matrix. You need the action of it. For example, check the "Using Function Handle Instead of Numeric Matrix" example at the link https://www.mathworks.com/help/matlab/ref/pcg.html . The main idea is that you create a function which takes a vector and computes A*x without constructing the matrix A. In your case, it could be a for loop.

The same idea applies to Jacobi, too. However, for ICHOL, you need the matrix explicitly and there is no way around it. Hence, my first point. Probably, you are not expected to use more than 500 points in each direction (or 250,000 points in total).

• I have tried to implement something like what they show in the matlab link you have provided me. I attached my attempt in the question as an edit, but when I do it that way it doesn't work for me, the error goes wrong, I don't know what I have wrong. And it is that they explicitly ask me to do it with 5,000 and 10,000 points per side ... – Haus Jun 26 at 6:19
• I tried over the weekend, 10000 points per side is do-able using matrices if you have around 16GB of RAM. But of course, you can't use a direct solver (except for forward and backward substitution for ichol). However, at 1000 points per side the error is already at 1e-8 range and using more than that will not (probably) improve that. Also, iterative solvers with the preconditioners you mentioned (as I would expect) are struggling to solve the linear system at those many points. I don't really understand the point of this exercise. – Abdullah Ali Sivas Jun 28 at 10:33
• If I have time sometime today, I will update my answer. – Abdullah Ali Sivas Jun 28 at 10:34