# Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

[u,ux,uy,X,Y]=laplarec(...
```function [u,ux,uy,X,Y]=laplarec(...
ubot,utop,ulft,urht,a,b,nx,ny,N)
%
% [u,ux,uy,X,Y]=laplarec(...
%               ubot,utop,ulft,urht,a,b,nx,ny,N)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This program evaluates a harmonic function and its
% first partial derivatives in a rectangular region.
% The method employs a Fourier series expansion.
% ubot   - defines the boundary values on the bottom
%          side. This can be an array in which
%          ubot(:,1) is x coordinates and ubot(:,2)
%          is function values. Values at intermediate
%          points are obtained by piecewise linear
%          interpolation. A character string giving
%          the name of a function can also be used.
%          Then the function is evualuated using 200
%          points along a side to convert ubot to an
%          array. Similar comments apply for utop,
%          ulft, and urht introduced below.
% utop   - boundary value definition on the top side
% ulft   - boundary value definition on the left side
% urht   - boundary value definition on the right side
% a,b    - rectangle dimensions in x and y directions
% nx,ny  - number of x and y values for which the
%          solution is evaluated
% N      - number of terms used in the Fourier series
% u      - function value for the solution
% ux,uy  - first partial derivatives of the solution
% X,Y    - coordinate point arrays where the solution
%          is evaluated
%
% User m functions used: datafunc ulinbc
%                        recseris ftsincof

disp(' ')
disp('SOLVING THE LAPLACE EQUATION IN A RECTANGLE')
disp(' ')

if nargin==0
disp(...
'Give the name of a function defining the data')
%datfun=input(...
%   '(try datafunc as an example): > ? ','s');
datfun='datafunc';
[ubot,utop,ulft,urht,a,b,nx,ny,N]=feval(datfun);
end

% Create a grid to evaluate the solution
x=linspace(0,a,nx); y=linspace(0,b,ny);
[X,Y]=meshgrid(x,y); d=(a+b)/1e6;
xd=linspace(0,a,201)'; yd=linspace(0,b,201)';

% Check whether boundary values are given using
% external functions. Convert these to arrays

if isstr(ubot)
ud=feval(ubot,xd); ubot=[xd,ud(:)];
end
if isstr(utop)
ud=feval(utop,xd); utop=[xd,ud(:)];
end
if isstr(ulft)
ud=feval(ulft,yd); ulft=[yd,ud(:)];
end
if isstr(urht)
ud=feval(urht,yd); urht=[yd,ud(:)];
end

% Determine function values at the corners
ub=interp1(ubot(:,1),ubot(:,2),[d,a-d]);
ut=interp1(utop(:,1),utop(:,2),[d,a-d]);
ul=interp1(ulft(:,1),ulft(:,2),[d,b-d]);
ur=interp1(urht(:,1),urht(:,2),[d,b-d]);
U=[ul(1)+ub(1),ub(2)+ur(1),ur(2)+ut(2),...
ut(1)+ul(2)]/2;

% Obtain a solution satisfying the corner
% values and varying linearly along the sides

[v,vx,vy]=ulinbc(U,a,b,X,Y);

% Reduce the corner values to zero to improve
% behavior of the Fourier series solution
% near the corners

f=inline('u0+(u1-u0)/L*x','x','u0','u1','L');
ubot(:,2)=ubot(:,2)-f(ubot(:,1),U(1),U(2),a);
utop(:,2)=utop(:,2)-f(utop(:,1),U(4),U(3),a);
ulft(:,2)=ulft(:,2)-f(ulft(:,1),U(1),U(4),b);
urht(:,2)=urht(:,2)-f(urht(:,1),U(2),U(3),b);

% Evaluate the series and combine results
% for the various component solutions

[ub,ubx,uby]=recseris(ubot,a,b,1,x,y,N);
[ut,utx,uty]=recseris(utop,a,b,2,x,y,N);
[ul,ulx,uly]=recseris(ulft,a,b,3,x,y,N);
[ur,urx,ury]=recseris(urht,a,b,4,x,y,N);
u=v+ub+ut+ul+ur; ux=vx+ubx+utx+ulx+urx;
uy=vy+uby+uty+uly+ury; close

% Show results graphically

surfc(X,Y,u), xlabel('x axis'), ylabel('y axis')
zlabel('U(X,Y)')
title('HARMONIC FUNCTION IN A RECTANGLE')
colormap(gray), brighten(.7)
shg, pause
print -deps laprecsr

colormap('default')
contour(X,Y,u,30); title('Contour Plot');
xlabel('x direction'); ylabel('y direction');
colorbar
colormap(gray), brighten(.7)
shg, pause
print -deps laprecnt

surf(X,Y,ux), xlabel('x axis'), ylabel('y axis')
zlabel('DU(X,Y)/DX')
title('DERIVATIVE OF U(X,Y) IN THE X DIRECTION')
view(225,20), colormap(gray), brighten(.7)
shg, pause
print -deps laprecdx

surf(X,Y,uy), xlabel('x axis'), ylabel('y axis')
zlabel('DU(X,Y)/DY')
title('DERIVATIVE OF U(X,Y) IN THE Y DIRECTION')
colormap(gray), brighten(.7), pause
print -deps laprecdy
shg

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

function [ubot,utop,ulft,urht,a,b,...
nx,ny,N]=datafunc
%
% [ubot,utop,ulft,urht,a,b,...
%         nx,ny,N]=datafunc
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This is a sample data case which can be
% modified to apply to other examples
%
% ubot, utop - vectors of function values on the
%              bottom and top sides
% ulft, urht - vectors of function values on the
%              right and left sides
% a, b       - rectangle dimensions along the
%              x and y axis
% nx, ny     - number of grid values for the x
%              and y directions
% N          - number of terms used in the
%              Fourier series solution

a=3; b=2; e=1e-5; N=100;
x=linspace(0,1,201)'; s=sin(pi*x);
c=cos(pi*x); ubot=[a*x,2-4*s];
utop=[a*x,interp1([0,1/3,2/3,1],...
[-2,2,2,-2],x)];
ulft=[b*x,2*c]; urht=ulft; nx=51; ny=31;

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

function [u,ux,uy]=ulinbc(U,a,b,X,Y)
%
% [u,ux,uy]=ulinbc(U,a,b,X,Y)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function determines a harmonic function
% varying linearly along the sides of a rectangle
% with specified corner values
%
% U     - corner values of the harmonic function
%         [U(1),...U(4)] <=> corner coordinates
%         (0,0), (0,a), (a,b), (0,b)
% a,b   - rectangle dimensions in the x and y
%         directions
% X,Y   - array coordinates where the solution
%         is evaluated
% u     - function values evaluated for X,Y
% ux,uy - first derivative components evaluated
%         for the X,Y arrays

c=[1,0,0,0;1,a,0,0;1,a,b,a*b;1,0,b,0;]\U(:);
u=c(1)+c(2)*X+c(3)*Y+c(4)*X.*Y;
ux=c(2)+c(4)*Y; uy=c(3)+c(4)*X;

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

function [u,ux,uy,X,Y]=recseris(udat,a,b,iside,x,y,N)
%
% [u,ux,uy,X,Y]=recseris(udat,a,b,iside,x,y,N)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes a function harmonic in
% a rectangle with general function values given
% on one side and zero function values on the
% other three sides.
% udat    - a data array to determine the function
%           values by piecewise linear interpolation
%           along the side having nonzero values.
%           udat(:,1) contains either x or y values
%           along a side, and udat(:,2) contains
%           corresponding function values
% a,b     - side lengths for the x and y directions
% iside   - an index indicating the side for which
%           function values are given.
%           [1,2,3,4]<=>[bottom,top,left,right]
% x,y       data vectors defining a grid
%           [X,Y]=meshgrid(x,y) on which the function
%           and its first partial derivatives are
%           computed
% N       - number of series terms used (up to 500)
% u,ux,uy - arrays of values of the harmonic function
%           and its first partial derivatives
% X,Y       arrays of coordinate values for which
%           function values were computed.

x=x(:)'; y=y(:); ny=length(y); N=min(N,500);
if iside<3, period=2*a; else, period=2*b; end
c=ftsincof(udat,period); n=1:N; c=c(n);
if iside<3     % top or bottom sides
npa=pi/a*n; c=c./(1-exp(-2*b*npa));
sx=sin(npa(:)*x); cx=cos(npa(:)*x);
if iside==1 % bottom side
dy=exp(-y*npa); ey=exp(-(2*b-y)*npa);
u=repmat(c,ny,1).*(dy-ey)*sx;
c=repmat(c.*npa,ny,1);
ux=c.*(dy-ey)*cx; uy=-c.*(dy+ey)*sx;
else        % top side
dy=exp((y-b)*npa); ey=exp(-(y+b)*npa);
u=repmat(c,ny,1).*(dy-ey)*sx;
c=repmat(c.*npa,ny,1);
ux=c.*(dy-ey)*cx; uy=c.*(dy+ey)*sx;
end
else           % left or right sides
npb=pi/b*n; c=c./(1-exp(-2*a*npb));
sy=sin(y*npb); cy=cos(y*npb);
if iside==3 % left side
dx=exp(-npb(:)*x);
ex=exp(-npb(:)*(2*a-x));
u=repmat(c,ny,1).*sy*(dx-ex);
c=repmat(c.*npb,ny,1);
ux=c.*sy*(-dx-ex); uy=c.*cy*(dx-ex);
else        % right side
dx=exp(-npb(:)*(a-x));
ex=exp(-npb(:)*(a+x));
u=repmat(c,ny,1).*sy*(dx-ex);
c=repmat(c.*npb,ny,1);
ux=c.*sy*(dx+ex); uy=c.*cy*(dx-ex);
end
end
[X,Y]=meshgrid(x,y);

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

function c=ftsincof(y,period)
%
% c=ftsincof(y,period)
% ~~~~~~~~~~~~~~~~~~~
% This function computes 500 Fourier sine
% coefficients for a piecewise linear
% function defined by a data array
% y      - an array defining the function
%          over half a period as
%          Y(x)=interp1(y(:,1),y(:,2),x)
% period - the period of the function
%
xft=linspace(0,period/2,513);
uft=interp1(y(:,1),y(:,2)/512,xft);
c=fft([uft,-uft(512:-1:2)]);
c=-imag(c(2:501));
```