Example 5.8

Iterative Solution of an Elliptic Equation

Contents

Setup

close all;

N = 20;  % points in the each direction

x0 = -1;
xN = 1;
y0 = -1;
yN = 1;

dx = (xN-x0)/N;
dy = (yN-y0)/N;

x = linspace(x0,xN,N+1);
y = linspace(y0,yN,N+1);

[X Y] = meshgrid(x,y);

Q = 2*(2-X.^2 - Y.^2);
Qint = Q(2:N,2:N);  % obtain interior
q = Qint(:);        % transform into a vector

Direct Solve

construct system. See Applied Numerical Linear Algebra by James Demmel footnote on page 275

TN = gallery('tridiag',N-1,-1,2,-1);
TNxN = kron(speye(N-1),TN) + kron(TN,eye(N-1));
tic
pd = (1/(dx^2)*TNxN)\q;
toc
Pd = zeros(N+1,N+1);
Pd(2:N,2:N) = reshape(pd,N-1,N-1);

Elapsed time is 0.000690 seconds.

True solution

Pd = (X.^2-1).*(Y.^2-1);

tol = 1e-4;

Point Jacobi

Pj = zeros(N+1,N+1);   % initial guess
Pnew = zeros(N+1,N+1); % extra storage
err = max(max(abs(Pd(2:N,2:N)-Pj(2:N,2:N))./Pd(2:N,2:N)));
iter = 0;

while ( err > tol )
    iter = iter + 1;
    for j = 2:N
        for i = 2:N
            Pnew(i,j) = (Pj(i+1,j)+Pj(i-1,j)+Pj(i,j+1)+Pj(i,j-1))/4 + dx^2/4*Q(i,j);
        end
    end
    Pj = Pnew;
    err = max(max(abs(Pd-Pj)))/max(max(abs(Pd)));
    Ej(iter) = err;
end

pj_iter = iter

pj_iter =

   749

Gauss-Seidel

Pgs = zeros(N+1,N+1);   % initial guess
err = max(max(abs(Pd(2:N,2:N)-Pgs(2:N,2:N))./Pd(2:N,2:N)));
iter = 0;

while ( err > tol )
    iter = iter + 1;
    for j = 2:N
        for i = 2:N
            Pgs(i,j) = (Pgs(i+1,j)+Pgs(i-1,j)+Pgs(i,j+1)+Pgs(i,j-1))/4 + dx^2/4*Q(i,j);
        end
    end
    err = max(max(abs(Pd-Pgs)))/max(max(abs(Pd)));
    Egs(iter) = err;
end

pgs_iter = iter

pgs_iter =

   375

SOR

Psor = zeros(N+1,N+1);   % initial guess
Pold = zeros(N+1,N+1);   % old solution
omega = 1.8; % SOR parameter
err = max(max(abs(Pd(2:N,2:N)-Psor(2:N,2:N))./Pd(2:N,2:N)));
iter = 0;

while ( err > tol )
    iter = iter + 1;
    for j = 2:N
        for i = 2:N
            Psor(i,j) = omega*((Psor(i+1,j)+Psor(i-1,j)+Psor(i,j+1)+Psor(i,j-1))/4 + dx^2/4*Q(i,j))+(1-omega)*Psor(i,j);
        end
    end
    err = max(max(abs(Pd-Psor)))/max(max(abs(Pd)));
    Esor(iter) = err;
end

psor_iter = iter

psor_iter =

    45

Published with MATLAB® 7.5