Search code examples
matlabmathequationpdemathematical-expressions

A problem in using Five-point difference format solve Poisson equation with Matlab


The problem set is this, a poisson equation. I want to solve it with Five-point difference format. But i have some problem in the code.

enter image description here

The code

function F = fivepointdiff(~,n)
h=1/n;
N=2*(n-1)*n+(3*n-1)*(n-1);
XY=zeros(2,n);
for i = 1:n
    for j=1:n-1
        XY(:,(n-1)*(i-1)+j)=[1+j*h;i*h];
    end
end
for i =1:n-1
    for j=1:3*n-1
        XY(:,n*(n-1)+(3*n-1)*(i-1)+j)=[j*h;1+i*h];
    end
end
for i=1:n
    for j=1:n-1
        XY(:,n*(n-1)+(3*n-1)*(n-1)+(n-1)*(i-1)+j)=[1+j*h;2*1+(i-1)*h];
    end
end
A=zeros(N,N);
for i=1:N
    for j=1:N
        if(i==j)
            A(i,j)=4;
        else if(((XY(1,i)-XY(1,j))^2+(XY(2,i)-XY(2,j))^2)<2*h*h)
                A(i,j)=-1;
            end
        end
    end
end
f=zeros(N,1);
for i =1:N
    f(i,1)=h*h;
end
U=bicg(A,f,0,1,100);
F=[XY;U'];
                

When I run it, There are some error

fivepointdiff(1, 25)

Warning: Input tol may not be achievable by BICG Try to use a bigger tolerance

In bicg (line 104)

In fivepointdiff (line 35)

Error using bicg (line 135)

Preconditioner must be a square matrix of size 2976 to match the problem size.

Error in fivepointdiff (line 35)

U=bicg(A,f,0,1,100);


Solution

  • You are using the function 'bicg' the wrong way.

    According to: https://www.mathworks.com/help/matlab/ref/bicg.html

    x = bicg(A,b,tol,maxit,M)
    
    • The third argument 'tol' is the tolerance and it can't be equal to 0 and this is the reason of the warning you get. (Warning: Input tol may not be achievable by BICG Try to use a bigger tolerance)
    • 'maxit' is the maximum number of iterations and it is not logical to set it to 1.
    • Finally, M should be a matrix with the same dimensions as A and not a number as you specified.Using this matrix can improve the numerical properties of the problem and the efficiency of the calculation and its use is optional.

    So, you need to modify the function arguments, maybe something like this:

    x = bicg(A,b,1e-6,1000)