Problem nr. 3

function U=lapsolve1(bot,right,top,left)
% Solve Laplace equation on a square.
% Parameters:
%
%    bot,right,top,left are vectors of boundary values,
%    each has length n=nr. of inner nodes/axis
% The solution is returned as an (n+2)x(n+2) matrix (the boundary
% values are appended)
n=length(bot);  % nr. of interior nodes
m=n+1;          % nr. of subintervals
% Assemly of A, use function lapm that has been written earlier.
A=lapm(m);       % Make sure to include this within Matlab path
% Assembly of RHS
B=[bot;[left(2:n-1)' zeros(n-2,n-2) right(2:n-1)'];top];
B(1,1)=bot(1)+left(1);B(1,n)=bot(n)+right(1);
B(n,1)=left(n)+top(1);B(n,n)=top(n)+right(n);
B=B';B=B(:);
% Now the system is ready, then just solve it.
uvek=A\B;
U=reshape(uvek,n,n)';
U=flipud(U);
% Append boundary values.

U=[fliplr(left)' U fliplr(right)'];

NW=0.5*(left(n)+top(1));NE=0.5*(top(n)+right(n));
SE=0.5*(bot(n)+right(1));SW=0.5*(bot(1)+left(1));

U=[NW top NE;
   U;
   SW bot SE];

% Remember to check everything went all right by
% del2

Code for lapm, used in lapsolve1

function A=lapm(m)
% Matrix for solving Laplace (Poisson) equation with 5 point. discr.
% Sparse structure.
% Arguments: m -- nr. of subintervals
%
n=m-1;  % 
N=n^2; % nr. of inner nodes = nr. of unknowns
e=ones(N,1);
d=~rem(1:N,n)==0;du=[1,d];du=du(1:N);
A=-spdiags([e,d',-4*e,du',e],[-n,-1,0,1,n],N,N);

First test

left=[20,40,60,80];top=[50,20,10,0];bot=zeros(1,4),right=left
U=lapsolve1(bot,right,top,left)
n=length(left);

for i=2:n+1; for j=2:n+1 ;
  d2U(i,j)= 4*U(i,j)-U(i-1,j)-U(i+1,j)-U(i,j-1)-U(i,j+1);
end; end
d2U
format long

d2U(2:n+1,2:n+1)
ans =

   1.0e-13 *

  -0.35527136788005  -0.07105427357601  -0.05329070518201                  0
                  0   0.07105427357601   0.00888178419700  -0.08881784197001
   0.02664535259100  -0.02664535259100  -0.02664535259100  -0.00888178419700
                  0   0.00888178419700                  0  -0.01776356839400

You can also use the Matlab function:

del2(U)

Finer grid, interpolate boundary data

 POLYFIT Fit polynomial to data.
    POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
    degree N that fits the data, P(X(I))~=Y(I), in a least-squares sense.
  POLYVAL Evaluate polynomial.
    Y = POLYVAL(P,X), when P is a vector of length N+1 whose elements
    are the coefficients of a polynomial, is the value of the
    polynomial evaluated at X.
 
        Y = P(1)*X^N + P(2)*X^(N-1) + ... + P(N)*X + P(N+1)
 
    If X is a matrix or vector, the polynomial is evaluated at all
    points in X. 

 xd=linspace(0,1,6);xd=xd(2:5);  % Given x (or y) datapoints, fixed.

 m=10;n=m-1;xi=linspace(0,1,m+1);
% To run this with various grid densities, you just need to change
% m=... and nothing else.

yd=left;
 c=polyfit(xd,yd,length(xd)-1)
 leftn=polyval(c,xi)
 plot(xd,left,'b*',xi,leftn,'ro')

yd=bot;
 c=polyfit(xd,yd,length(xd)-1)
 botn=polyval(c,xi);
 plot(xd,bot,'b*',xi,botn,'ro') 

yd=right;
 c=polyfit(xd,yd,length(xd)-1)
 rightn=polyval(c,xi);
 plot(xd,right,'b*',xi,rightn,'ro')

yd=top;
 c=polyfit(xd,yd,length(xd)-1)
 topn=polyval(c,xi);
 plot(xd,top,'b*',xi,topn,'ro')

Un=lapsolve1(botn,rightn,topn,leftn);

 surf(Un)
 figure;
 surf(U)

Next use 2-dim interpolation to compare error behaviour