No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

Howard Wilson

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

[u,ux,uy,X,Y]=recseris(udat,a,b,iside,x,y,N)
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);

Contact us