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.

vvheun3D(varargin)
```function [w1, w2, w3] = vvheun3D(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 1by3) - upper bound of spatial (x) variable
%       (Default L = [1,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 1by3) - number of discrete spatial intervals
%       (Default m = [100,100,100])
%   n (scalar) - number of discrete time intervals
%       (Default n = 100)
%
%   w (m+1 x 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);
L_z = L(3);
m_x = m(1);
m_y = m(2);
m_z = m(3);

% initialize h, k, lambda, and w
h_x = L_x/m_x;
h_y = L_y/m_y;
h_z = L_z/m_z;
k = T./n;
w1 = zeros(m_y+1,m_x+1,m_z+1,n);
w2 = zeros(m_y+1,m_x+1,m_z+1,n);
w3 = zeros(m_y+1,m_x+1,m_z+1,n);

% initialize w1 (to sin(pix1/L1)*sin(pix2/L2)*sin(pix3/L3))
for i=2:m_x
for j=2:m_y
for k=2:m_z

w1(j,i,k,1) = sin(pi.*h_x.*(i-1)./L_x).*sin(pi.*h_y.*(j-1)./L_y)...
.*sin(pi.*h_z.*(k-1)./L_z);

end
end
end

% initialize w2 (to 4*sin(pix1/L1)*sin(pix2/L2)*sin(pix3/L3))
for i=2:m_x
for j=2:m_y
for k=2:m_z

w2(j,i,k,1) = 4*sin(pi.*h_x.*(i-1)./L_x).*sin(pi.*h_y.*(j-1)./L_y)...
.*sin(pi.*h_z.*(k-1)./L_z);

end
end
end

% initialize w3 (to sin(pix1/L1)*sin(pix2/L2)*sin(pix3/L3))
for i=2:m_x
for j=2:m_y
for k=2:m_z

w3(j,i,k,1) = 2*sin(pi.*h_x.*(i-1)./L_x).*sin(pi.*h_y.*(j-1)./L_y)...
.*sin(pi.*h_z.*(k-1)./L_z);

end
end
end

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

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

g3 = zeros(m_y+1,m_x+1,m_z+1);
p3 = zeros(m_y+1,m_x+1,m_z+1);
o3 = zeros(m_y+1,m_x+1,m_z+1);

for t=1:(n-1)
%Go back and check that these coefficients are correct.
f1 = (alpha^2)*laplacian(w1(:,:,:,t),h_y,h_x,h_z);
g1 = w1(:,:,:,t) + ((2*k)/3)*f1;
p1 = (alpha^2)*laplacian(g1,h_y,h_x,h_z);
o1 = p1 + ((2*k)/3)*(alpha^2)*laplacian(p1,h_y,h_x,h_z);

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

f3 = (alpha^2)*laplacian(w3(:,:,:,t),h_y,h_x,h_z);
g3 = w3(:,:,:,t) + ((2*k)/3)*f3;
p3 = (alpha^2)*laplacian(g3,h_y,h_x,h_z);
o3 = p3 + ((2*k)/3)*(alpha^2)*laplacian(p3,h_y,h_x,h_z);

w1(:,:,:,t+1) = w1(:,:,:,t) + (k/4)*f1 + ((3*k)/4)*o1;
w2(:,:,:,t+1) = w2(:,:,:,t) + (k/4)*f2 + ((3*k)/4)*o2;
w3(:,:,:,t+1) = w3(:,:,:,t) + (k/4)*f3 + ((3*k)/4)*o3;
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,100];
end

% check m has valid value
if (numel(m)~=3) || 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,1];
end

% check L has valid value
if (numel(L)~=3) || any(L<=0)
error('L must be a positive vector with 3 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.');
end

%another subfunction
function y = laplacian(x,h_y,h_x,h_z)
[m,n,p] = size(x);
y = zeros(m,n,p);
for i = 2:n-1
for j = 2:m-1
for k = 2:p-1

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

end
end
end```