Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

laplarec
function laplarec
% Example:  laplarec
% ~~~~~~~~~~~~~~~~~~
% This program uses Fourier series methods to
% solve Laplace's equation inside a rectangle.
% Boundary conditions are defined by piecewise
% linear interpolation of boundary data. The
% program can easily be changed to deal with
% problems where the boundary conditions are
% expressed by analytically defined functions.
% Surface and contour plots of the solution
% are also provided.
%
% User m functions required:
%    setup, ulbc, sinfft, laprec, fbot, gtop,
%    plft, qrht, lintrp

clear
global a_ b_ xbot_ ubot_ xtop_ utop_
global ylft_ ulft_ yrht_ urht_
global u1_ u2_

a_=3; b_=2; % Rectangle side lengths

a=a_; b=b_;
xd=linspace(0,a,101)'; yd=linspace(0,b,101)';
fb=[xd,1+2*(xd/a)+sin(4*pi/a*xd)/2]; % cos(5*pi/a*xd)];
ft=[a*[0;.3;.4;.6;.7;1],[2;2;3;3;2;2.5]];
fl=[[0;b/2;b],[1;2;2]]; fr=[[0;b/5;b],[3;3;2.5]];
xtop_=ft(:,1); utop_=ft(:,2); ntop=length(xtop_);
xbot_=fb(:,1); ubot_=fb(:,2); nbot=length(xbot_);
ylft_=fl(:,1); ulft_=fl(:,2); nlft=length(ylft_);
yrht_=fr(:,1); urht_=fr(:,2); nrht=length(yrht_);

% Data for piecewise linear boundary conditions
while 0
xtop_=[0;    3];    utop_=  [0;    2];
ntop=length(xtop_);
xbot_=[0;1;2;3];    ubot_=[2;-1;-1;2];
nbot=length(xbot_);
ylft_=[0;0.6;0.6;2]; ulft_=[2;2;4;-2];
nlft=length(ylft_);
yrht_=[0;1;1.1;2];  urht_=[2;4;-2; 0];
nrht=length(ylft_);
end

% Create the constants needed in function ulbc
cc=setup;

% Adjust boundary data to give zero
% end conditions
for k=1:nbot
ubot_(k)=ubot_(k)-ulbc(cc,xbot_(k),0);end
for k=1:ntop
utop_(k)=utop_(k)-ulbc(cc,xtop_(k),b_);end
for k=1:nlft
ulft_(k)=ulft_(k)-ulbc(cc,0,ylft_(k));end
for k=1:nrht
urht_(k)=urht_(k)-ulbc(cc,a_,yrht_(k));end

% Generate Fourier coefficients for the
% modified boundary conditions
sbot=sinfft('fbot',a_,9); stop=sinfft('gtop',a_,9);
slft=sinfft('plft',b_,9); srht=sinfft('qrht',b_,9);

% Generate a grid of interior points
% for solution evaluation
nin=51; ntrms=50;
xin=linspace(0,a_,nin); yin=linspace(0,b_,nin);

% Evaluate the solution having zero
% corner values
uin=laprec(...
sbot,stop,slft,srht,a_,b_,xin,yin,ntrms);

% Correct results for nonzero corner values

uin=uin+ulbc(cc,xin,yin); uin=flipud(uin');

% Display surfaces showing function
% values on the grid

clf; surfc(xin,yin,flipud(uin));
xlabel('x axis'); ylabel('y axis');
zlabel('function value');
title('Surface Plot'); axis equal, figure(gcf);
dumy=input('Press [Enter] to continue','s');
print -deps lapsrfac

clf; contour(xin,yin,flipud(uin),30);
title('Contour Plot');
xlabel('x direction'); ylabel('y direction');
figure(gcf), print -deps contur
disp('All Done');

%=============================================

function cc=setup
%
% cc=setup
% ~~~~~~~~
% cc - coefficients used to define the
%      particular solution to satisfy
%      corner conditions
%
% User m functions called:
%    fbot, qrht, plft, gtop
%----------------------------------------------

global  a_ b_

s=(a_+b_)/1e10;
c1=(fbot(s)+plft(s))/2;
c2=(fbot(a_-s)+qrht(s))/2;
c3=(plft(b_-s)+gtop(s))/2;
c4=(gtop(a_-s)+qrht(b_-s))/2;
mat=[[1,0,0,0];[1,a_,0,0];
[1,0,b_,0];[1,a_,b_,a_*b_]];
vec=[c1;c2;c3;c4]; cc=mat\vec;

%=============================================

function u=ulbc(c,x,y)
%
% u=ulbc(c,x,y)
% ~~~~~~~~~~~~~
% This function determines a harmonic function
% satisfying boundary conditions which vary
% linearly on the sides of a rectangle.
%
% User m functions called:  none
%----------------------------------------------

x=x(:); y=y(:)'; nx=length(x); ny=length(y);
x=x*ones(1,ny); y=ones(nx,1)*y;
u=c(1)*ones(nx,ny)+c(2)*x+c(3)*y+c(4)*x.*y;
function sincof = sinfft(fun,hafper,powr2)
%
% sincof = sinfft(fun,hafper,powr2)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function determines coefficients in
% the Fourier sine series of a general real
% valued function.
%
%  hafper - the half  period over which the
%           expansion applies
%  fun    - the real valued function being
%           expanded. This function must be
%           defined for vector arguments with
%           components between zero
%           and hafper.
%  powr2  - the power of 2 determining the
%           number of function values used
%           in the FFT (number of
%           values = 2^powr2).  When powr2
%           is 9, then 512 function values
%           are used and 255 Fourier
%           coefficients are computed.
%
% User m functions called:  none
%----------------------------------------------

n=2^powr2; period=2*hafper; n2=n/2;
x=(period/n)*(0:n2);
fval=feval(fun,x);
fval=fval(:); fval=[fval;-fval(n2:-1:2)];
foucof=fft(fval);
sincof=-(2/n)*imag(foucof(2:n2));

%=============================================

function u=laprec(f,g,p,q,a,b,x,y,ntrms)
%
% u=laprec(f,g,p,q,x,y,ntrms)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function sums the series which solves
% Laplace's equation in a rectangle.
%
%  f,g,p,q - Fourier sine series coefficients
%            for the boundary conditions on
%            the bottom, top, left, and
%            right sides
%   a,b    - the horizontal and vertical
%            side lengths
%   x,y    - vectors containing coordinates
%            of points defining a rectangular
%            grid on which the solution is
%            to be evaluated.
%   ntrms  - number of series terms used (not
%            exceeding length(f));
%
% User m functions called:  none
%----------------------------------------------

nt=length(f);
if nargin==8, ntrms=nt; end;
ntrms=min(nt,ntrms);
x=x(:); y=y(:)';
n=1:ntrms; nx=length(x); ny=length(y);
if nt>ntrms
f=f(n); g=g(n); p=p(n); q=q(n);
end
a2=2*a; b2=2*b; na=(pi/a)*n; nb=(pi/b)*n;
denomx=1-exp(-b2*na(:));
f=f(:)./denomx; g=g(:)./denomx;
denomy=1-exp(-a2*nb(:)');
p=p(:)'./denomy; q=q(:)'./denomy;

u1_=(f*ones(1,ny)).* ...
(exp(-na'*y)-exp(-na'*(b2-y)));
u2_=(g*ones(1,ny)).* ...
(exp(-na'*(b-y))-exp(-na'*(b+y)));
u=sin(x*na)*( u1_+u2_ );
u3_=(exp(-x*nb)-exp(-(a2-x)*nb)).* ...
(ones(nx,1)*p);
u4_=(exp(-(a-x)*nb)-exp(-(a+x)*nb)).* ...
(ones(nx,1)*q);
u=u+(u3_+u4_)*sin(nb'*y);

%=============================================

function ubot=fbot(x)
%
% ubot=fbot(x)
% ~~~~~~~~~~~~
% x    - vector argument
% ubot - function value on bottom side
%
% User m functions called:  lintrp
%----------------------------------------------

global xbot_ ubot_

ubot=lintrp(xbot_,ubot_,x);

%=============================================

function utop=gtop(x)
%
% utop=gbot(x)
% ~~~~~~~~~~~~
% x    - vector argument
% gtop - function value on top side
%
% User m functions called:  lintrp
%----------------------------------------------

global xtop_ utop_

utop=lintrp(xtop_,utop_,x);

%=============================================

function ulft=plft(y)
%
% ulft=plft(y)
% ~~~~~~~~~~~~
% y    - vector argument
% ulft - function value on left side
%
% User m functions called:  lintrp
%----------------------------------------------

global ylft_ ulft_

ulft=lintrp(ylft_,ulft_,y);

%=============================================

function urht=qrht(y)
%
% urht=qrht(y)
% ~~~~~~~~~~~~
% y    - vector argument
% urht - function value on right side
%
% User m functions called:  lintrp
%----------------------------------------------

global yrht_ urht_

urht=lintrp(yrht_,urht_,y);