Code covered by the BSD License

Highlights fromMatlab code for my Graduate Thesis

from Matlab code for my Graduate Thesis by Troy
Numerically solves the diffusion equations as it pertains to medical imaging.

vvheun2Dbending(varargin)
```function [w1, w2] = vvheun2Dbending(varargin)
% crankNicolson: uses Crank-Nicolson algorithm to approximate the solution
%   to the parabolic PDE:
%       u_{t}(x,t) - alpha^2 u_{xx}(x,t) = 0, 0<x<L, 0<t<T
%   subject to the boundary conditions
%       u(0,t) = u(L,t) = 0, 0<t<T
%   and the initial conditions
%       u(x,0) = f(x), 0<=x<=L
%
% arguments:
%   L (vector 1by2) - upper bound of spatial (x) variable
%       (Default L = [1,1])
%   T (scalar) - upper bound of time (t) variable
%       (Default T = 1)
%   alpha (scalar) - square root of coefficient of u_{xx} term
%       (Default alpha = 1)
%   m (vector 1by2) - number of discrete spatial intervals
%       (Default m = [100,100])
%   n (scalar) - number of discrete time intervals
%       (Default n = 100)
%
%   w (m+1 x m+1 x n) - approximation to u(x,t) at discrete space/time positions
%

% author: Troy J. Winkstern
% email: tjw8191@rit.edu
% date: 30 Jan 2011

% parse input arguments
[L,T,alpha,m,n] = parseInputs(varargin{:});
L_x = L(1);
L_y = L(2);
m_x = m(1);
m_y = m(2);

% initialize h, k, lambda, and w
h_x = L_x/m_x;
h_y = L_y/m_y;
k = T./n;
lambda_x = (alpha.^2).*k./(h_x.^2);
lambda_y = (alpha.^2).*k./(h_y.^2);
w1 = zeros(m_y+1,m_x+1,n);
w2 = zeros(m_y+1,m_x+1,n);

% initialize w1
for i=1:m_x+1
for j=2:m_y
w1(j,i,1) = cos(pi.*h_x.*(i-1)./L_x).*sin(pi.*h_y.*(j-1)./L_y);
end
end

% initialize w2
for i=2:m_x
for j=1:m_y+1
w2(j,i,1) = sin(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y);
end
end

g1 = zeros(m_y+1,m_x+1);
p1 = zeros(m_y+1,m_x+1);
o1 = zeros(m_y+1,m_x+1);

g2 = zeros(m_y+1,m_x+1);
p2 = zeros(m_y+1,m_x+1);
o2 = zeros(m_y+1,m_x+1);

for t=1:(n-1)
f1 = laplacian1(w1(:,:,t),h_y,h_x);
g1 = w1(:,:,t) + ((2*k)/3)*f1;
p1 = laplacian1(g1,h_y,h_x);
o1 = p1 + ((2*k)/3)*laplacian1(p1,h_y,h_x);

f2 = laplacian2(w2(:,:,t),h_y,h_x);
g2 = w2(:,:,t) + ((2*k)/3)*f2;
p2 = laplacian2(g2,h_y,h_x);
o2 = p2 + ((2*k)/3)*laplacian2(p2,h_y,h_x);

w1(:,:,t+1) = w1(:,:,t) + (k/4)*f1 + ((3*k)/4)*o1;
w2(:,:,t+1) = w2(:,:,t) + (k/4)*f2 + ((3*k)/4)*o2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction parseInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,T,alpha,m,n] = parseInputs(varargin)

% check number of input arguments
nargs = length(varargin);
error(nargchk(0,5,nargs));

% get/set n
if nargs<5
n = [];
else
n = varargin{5};
end
if isempty(n)
n = 100;
end

% check n has valid value
if (numel(n)>1) || (n<1) || ~isequal(round(n),n)
error('n must be a positive integer.');
end

% get/set m
if nargs<4
m = [];
else
m = varargin{4};
end
if isempty(m)
m = [100,100];
end

% check m has valid value
if (numel(m)~=2) || any(m<1) || any(~isequal(round(m),m))
error('m must be a vector of positive integers.');
end

% get/set alpha
if nargs<3
alpha = [];
else
alpha = varargin{3};
end
if isempty(alpha)
alpha = 1;
end

% check alpha has valid value
if numel(alpha)>1
error('alpha must be a real scalar.');
end

% get/set T
if nargs<2
T = [];
else
T = varargin{2};
end
if isempty(T)
T = 1;
end

% check T has valid value
if (numel(T)>1) || (T<=0)
error('T must be a positive real scalar.');
end

% get/set L
if nargs<1
L = [];
else
L = varargin{1};
end
if isempty(L)
L = [1,1];
end

% check L has valid value
if (numel(L)~=2) || any(L<=0)
error('L must be a positive vector with 2 elements.');
end

% warning if L not integer
if any(~isequal(round(L),L))
warning('L should be an vector of integers so that the inital condition takes on the value of 0 for t = 0 and t = T.');
endarning('L should be an integer so that the inital condition takes on the value of 0 for t = 0 and t = T.');
end

%another subfunction
function y = laplacian1(x,h_y,h_x)
[m,n] = size(x);
y = zeros(m,n);

for j = 2:m-1
y(j,1) = (2*x(j,2)-2*x(j,1))/(h_x.^2)+(x(j+1,1)-2*x(j,1)+x(j-1,1))/(h_y.^2);
end

for i = 2:n-1
for j = 2:m-1

y(j,i) = (x(j-1,i)+x(j+1,i)-2*x(j,i))/(h_y.^2)+(x(j,i-1)+x(j,i+1)-2*x(j,i))/(h_x.^2);

end
end

for j = 2:m-1
y(j,n) = (2*x(j,n-1)-2*x(j,n))/(h_x.^2)+(x(j+1,n)-2*x(j,n)+x(j-1,n))/(h_y.^2);
end

%yet another subfunction
function y = laplacian2(x,h_y,h_x)
[m,n] = size(x);
y = zeros(m,n);

for i = 2:n-1
y(1,i) = (x(1,i+1)-2*x(1,i)+x(1,i-1))/(h_x.^2)+(2*x(2,i)-2*x(1,i))/(h_y.^2);
end

for i = 2:n-1
for j = 2:m-1

y(j,i) = (x(j-1,i)+x(j+1,i)-2*x(j,i))/(h_y.^2)+(x(j,i-1)+x(j,i+1)-2*x(j,i))/(h_x.^2);

end
end

for i = 2:n-1
y(m,i) = (x(m,i+1)-2*x(m,i)+x(m,i-1))/(h_x.^2)+(2*x(m-1,i)-2*x(m,i))/(h_y.^2);
end
```