Code covered by the BSD License  

Highlights from
Matlab code for my Graduate Thesis

Matlab code for my Graduate Thesis

by

 

Numerically solves the diffusion equations as it pertains to medical imaging.

crankNicolson(varargin)
function w = crankNicolson(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 (scalar) - upper bound of spatial (x) variable
%       (Default L = 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 (scalar) - number of discrete spatial intervals
%       (Default m = 100)
%   n (scalar) - number of discrete time intervals
%       (Default n = 100)
%
%   w (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{:});

% initialize h, k, lambda, and w
h = L./m;
k = T./n;
lambda = (alpha.^2).*k./(h.^2);
w = zeros(m+1,n);

% initialize first column of w
w(2:m,1) = sin(pi.*h.*(1:(m-1)).');

% construct matrices A and B based on lambda
Adiags = repmat([-lambda./2 1+lambda -lambda./2],m-1,1);
A = spdiags(Adiags,[-1 0 1],m-1,m-1);
Bdiags = repmat([lambda./2 1-lambda lambda./2],m-1,1);
B = spdiags(Bdiags,[-1 0 1],m-1,m-1);

% loop through columns of w, solving linear system:
% A w(:,i+1) = B w(:,i)
for i=1:(n-1)
    w(2:m,i+1) = A\(B*w(2:m,i));
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;
end

% check m has valid value
if (numel(m)>1) || (m<1) || ~isequal(round(m),m)
    error('m must be a positive integer.');
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;
end

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

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

Contact us