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 = zeros(N,5); Ad(:,1) = -k(X(:),Y(:)+dx/2); Ad(:,2) = -k(X(:)+dx*0.5,Y(:)); Ad(:,4) = -k(X(:)-dx*0.5,Y(:)); Ad(:,5) = -k( X(:), Y(:)-dx/2 ); Ad(:,3) = -sum( Ad(:,[1,2,4,5]), 2 ); idx = find(isDirichlet(:)); L=Ad(:,1);L(idx)=;B=Ad(:,2);B(idx)=; C=Ad(:,3);C(idx)=;U=Ad(:,4); U(idx)=;UP=Ad(:,5);UP(idx)=; 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