Solve boundary value problem by command lines

Hongxue Cai (view profile)

30 Jan 2006 (Updated )

solve PDEs in continous domain using command lines

solveBV
```function solveBV
% Codes showing how to solve a boundary value problem by command lines
%
%   1.  Construct geometry by a geometrical matrix
%   2.  Define boundary conditions (BCs) by a BC matrix
%   3.  Match coefficients between the MATLAB elliptic solver and the PDE to
%       be solved
%   4.  Discretize PDE and BCs in the geometrical domain
%   5.  Solve the BC problem and do post-processing
%
% These code need MATLAB PDE Toolbox
%
%       By Hongxue Cai  (h-cai@northwestern.edu)
%

% ===============  Problem to be solved   ===============>
%
% An incompressible fluid is in the domain deifned by the geometrical
% matrix gGeo below. This 2D domain is a cross section of a fluid chamber
% (scala tympani)
% in the cochlea, a spiral fluid-solid interacted structure in the inner ear.
% A deformable boundary (geometrical element #2 defined by the column# 2 in gGeo)
% represents a vibrating membrane (the basilar membrane). We want to find
% the distribution of the fluid pressure when a motion of the memebrane is
% given
%
% Equation to be solved:     L(P) = k^2* P          (1)
% where     P: fluid pressure
%           L: represent the Laplace operator
%           k: the wave number of the traveling wave on the memebrane
%
% Note: This equation is derived from the Navier-Stokes equation and the
% incomressibility of the fluid. In fact, the P is the dominant term of the
% WKB method applied to the pressure (see references)
%
% Bounary conditions:
%       dP/dn = r*w^2*U     on deformable boundary      (2)
%       dP/dn = 0           on rigid walls              (3)
% where:    r: the intensity of the fluid
%           n: the normal direction on a boundary point
%           w = 2*Pi*freq, with freq, the vibrating frequency
%           U:  normal displacement component in the direction of n
%
% MATLAB solver used: -grad(c*grad(P)) + a*P = f, with the following genral
% Neuman condtion:      n.(c*grad(P)) + q*P = g
% where c, a, f, q, g are coefficients to be defined by matching the PDE
% equatioin and BCs with the MATLAB solver, and we have
%       c = -1
%       a = k^2
%       f = 0
%       q = 0
%       g = 0 or -r*w^2*U
%
%References:
% 1. Cai, H., Manoussaki, D. and Chadwick, R.S., Effects of coiling on the micromechanics
%       of the mammalian cochlea, J. R. Soc. Interface, 2(4):341-348, 2005.
% 2. Cai, H., Shoelson, B and Chadwick R.S., Evidence of tectorial membrane radial motion
%       in a propagating mode of a complex cochlear model,
%       Proc. Natl. Acad. Sci. USA, 101(16):6243-6248, 2004.
% 3. Cai, H. and Chadwick, R.S., Radial structure of traveling waves in the inner ear,
%       SIAM J. Appl. Math., 63(4):1105-1120, 2003.
%
% < ======================================================

% Geometrical metrix        -------------------------------------------->
BM1 = -0.0079;
BM2 = 0.0171;
R0 = 0.022;
xc = 0;
yc = 0;
yBM = -R0/2.5;
xR = -sqrt(R0*R0 - yBM*yBM);
xL = -xR;
gGeo = [
2       2	    2	    1
xR      BM1	    BM2	    xR
BM1	    BM2	    xL      xL
yBM     yBM     yBM     yBM
yBM     yBM     yBM     yBM
0	    0	    0	    1
1	    1	    1	    0
0       0       0       xc
0       0       0       yc
0       0       0       R0
];
%geometry matrix, in which, each column defines a geometrical element that
%can be either a line or a circular arc. (Note that any complex geometry can
%be composed by line and circular arts.) In each column, the number in row#1
%indicates the type of the geometrica element, 2 for line element, 1 for
%circular arts. Then, row#2 is the x-coordinate of the starting point;
%row#3, that of the ending point; Row #4 and #5 are, respectively the
%y-coordinates of the the starting and ending point. The number in the 6th
%anf 7th
%rows indicate the location of geometric domain defined, (0,1) on the left side of this
%geometric element (When you walk from the begining point to the ending point along the element,
%the domain is on your left), (1,0) on th right. The rows 8, 9, 10 are zero for line,
%and repesent the x, y-coordinates and the radius of the center point for a
%circular arc.
%
%The values that define the geometry come from an application which models
%the baislar
%membrane in the inner ear. These values should be different for different
%applications, therefore one should not pay much attention to the
%individual numbers themselves.

% plotting geometry
figure
pdegplot(gGeo)      %MATLAB PDE Toolbox needed

%meshing
[p,e,t]=initmesh(gGeo);             %initial mesh
[p,e,t]=refinemesh(gGeo,p,e,t);     % refine mesh
[p,e,t]=refinemesh(gGeo,p,e,t);
%meshing data are stored in p, e, t

% display mesh
figure
pdemesh(p,e,t);                  %PDE Toolbox needed

% < -------------------------------------------

% -----------------   Define constant or known variable -------------->
freq = 100;     %frequency, Hz
w = 2*pi*freq;
r = 1;          %fluid intensity, g/cm^3
rw2 = r*w*w;
%
k = 6 - 0.03i;   %complex wavenumber
%Note that the geometrical domain (the cross section) is in the x-y plane, and the wave
%propogates in the z- direction with a wavenumber k. The geometry of the
%cochlea chages slowly in the z-direction, and the WKB method allows us to
%solve approximately a 3D flow in a 2D doamin (See references)
%
BMa = pi*(BM1+BM2)/2/(BM1-BM2);	%no initial mode change
BMb = -pi/(BM1-BM2);		        %no initial mode change
Am = 1e-9;	%Amplitude, 1 nm. Initially, the defomable element (the basilar memebrane) is given a motion
%           y = Am*cos(BMa + BMb*x), a clamped beam with maximum
%           displacement at the middle. When x = BM1, BM2, y = 0 (BMa +
%           BMb*x = -pi/2, pi/2)
% < --------------------------------------------------------------------

% define PDE coefficients    --------------------->
c = -1;
a = k^2;
f = 0;
% < -------------------------------------------

% Define boudary condition  ------------------------------->
%
% In this Ex, the geometrical element #2 (define by the column#2 in the
% matrix gGeo, subject to a non-zero Neumann BC; all other boundaries (#1,
% 3, and 4) are rigid walls, therefore a zero Neumann BC applied to.
%
% How many geometrical elements
n_Geo = size(gGeo,2);

n_mb = 2;	%No. of deformable boundary in gGeo

q = 48;      % q = 0

% g 4 deformable boundary
clear g
temp = [mat2str(-rw2*Am),'*cos(x*',mat2str(BMb),'+',mat2str(BMa),')'];
g(1:length(temp),1)=(temp*1)';
% to generate special codes in the BC matrix

% boundary condition matrix
bc = bcNeumann(n_mb, q, g, n_Geo);

%solve the PDE
disp('Elliptic problem (fluid in the lower chamber....>')
press = assempde(bc, p, e, t, c, a, f);

% post-processing ...........
%plotting the pressure distribution
figure
pdeplot(p, e, t, 'xydata', real(press))

%  ============ --------- end of sloveBV   -------------========
%
%
%sub-function
function bc = bcNeumann(n_mb, q, g, nCol)
%	-------- Neumann boundary condition matrix (scalar case) -------
%	n_mb:	Non-zero Neumann boundary segments (moving segments) index
%	nCol:	Total number of  boundary segments, nCol-length(n_mb)=
%           zero Neumann boundary segments number
%
%	------- n.(c grad u)+ q u = g  -----------------
%
Lmb = length(n_mb); % how many deformable boundaries
%
Lq = length(q);
Lg = length(g);
%
% number of rows in boundary condition matrix bc
row_b = max(Lq + Lg);
row_b = row_b + 4;
%
% set default values to 0
bc = zeros(row_b, nCol);
%
bc(1,:)=1;	%dimension of the system N = 1 (scalar case)
bc(2,:)=0;	%Number of Dirichlet boundary conditions M = 0
bc(3,:)=1;	%Lq = 1 for zero Neumann boundary segments (q = 0)
bc(4,:)=1;	%Lg = 1 for zero Neumann boundary segments (g = 0)
%
%Lq for non-zero Neumann boundary segments
for ij=1:Lmb
bc(3,n_mb(ij)) = Lq(ij);
end

%Lg for non-zero Neumann boundary segments
for ij=1:Lmb
bc(4,n_mb(ij)) = Lg(ij);
end

bc(5:6,:)=48;	 %q = g = 0 (48) for zero Neumann boundary segments

%q for non-zero Neumann boundary segments
nRow(1:Lmb) = 4;
for ij=1:Lmb
bc(nRow(ij)+1:nRow(ij)+Lq(ij),n_mb(ij)) = q(1:Lq(ij),ij);
end
nRow=nRow+Lq;

%g for non-zero Neumann boundary segments
for ij=1:Lmb
bc(nRow(ij) + 1:nRow(ij) + Lg(ij), n_mb(ij)) = g(1:Lg(ij),ij);
end
nRow = nRow + Lg;
row_b = length(bc);
%
return

```