function [w1, w2, w3] = vvmodeuler3Dneumann(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 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;
lambda_x = (alpha.^2).*k./(h_x.^2);
lambda_y = (alpha.^2).*k./(h_y.^2);
lambda_z = (alpha.^2).*k./(h_z.^2);
w1 = zeros(2*m_y+1,2*m_x+1,2*m_z+1,n);
w2 = zeros(2*m_y+1,2*m_x+1,2*m_z+1,n);
w3 = zeros(2*m_y+1,2*m_x+1,2*m_z+1,n);
% initialize w1
for i=1:m_x+1
for j=1:m_y+1
for k=1:m_z+1
w1(j,i,k,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y)...
.*cos(pi.*h_z.*(k-1)./L_z);
end
end
end
% initialize w2
for i=1:m_x+1
for j=1:m_y+1
for k=1:m_z+1
w1(j,i,k,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y)...
.*cos(pi.*h_z.*(k-1)./L_z);
end
end
end
% initialize w3
for i=1:m_x+1
for j=1:m_y+1
for k=1:m_z+1
w1(j,i,k,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y)...
.*cos(pi.*h_z.*(k-1)./L_z);
end
end
end
for t=1:(n-1)
f1 = laplacian(w1(:,:,:,t),h_y,h_x,h_z);
y1 = w1(:,:,:,t) + (k/2)*f1;
f2 = laplacian(w2(:,:,:,t),h_y,h_x,h_z);
y2 = w2(:,:,:,t) + (k/2)*f2;
f3 = laplacian(w3(:,:,:,t),h_y,h_x,h_z);
y3 = w3(:,:,:,t) + (k/2)*f3;
% went out to length 2m
X1 = [sqrt(1./(2*m_z)).*real(fft(y1, 2*m_y,1)):zeros(1,201,1)];
X2 = [sqrt(1./(2*m_z)).*real(fft(y2, 2*m_y,1)):zeros(1,201,1)];
X3 = [sqrt(1./(2*m_z)).*real(fft(y3, 2*m_y,1)):zeros(1,201,1)];
Y1 = [sqrt(1./(2*m_z)).*real(fft(X1, 2*m_x,2)):zeros(201,1,1)];
Y2 = [sqrt(1./(2*m_z)).*real(fft(X2, 2*m_x,2)):zeros(201,1,1)];
Y3 = [sqrt(1./(2*m_z)).*real(fft(X3, 2*m_x,2)):zeros(201,1,1)];
V1 = [sqrt(1./(2*m_z)).*real(fft(Y1, 2*m_z,3)):zeros(1,1,201)];
V2 = [sqrt(1./(2*m_z)).*real(fft(Y2, 2*m_z,3)):zeros(1,1,201)];
V3 = [sqrt(1./(2*m_z)).*real(fft(Y3, 2*m_z,3)):zeros(1,1,201)];
% angles divided by m
a_x = repmat((pi./m_x).*(0:m_x),[m_y+1,1,m_z+1]);
a_y = repmat((pi./m_y).*(0:m_y).',[1,m_x+1,m_z+1]);
a_z = repmat(reshape((pi./m_z).*(0:m_z),[1,1,m_z+1]),[m_y+1,m_x+1,1]);
% changed c to accomodate for 3-D
c1 = 1 + ((lambda_x)*(lambda_y)*(lambda_z)*(h_x)^2*(h_y)^2*(h_z)^2) - ...
(1/2)*((lambda_x)*cos(a_x)*(h_x)^2).*((lambda_y)*cos(a_y)*(h_y)^2) ...
.*((lambda_z)*cos(a_z)*(h_z)^2) - (1/2)*(((lambda_x)*cos(a_x)*(h_x)^2).^2)...
.*(((lambda_y)*cos(a_y)*(h_y)^2).^2).*(((lambda_z)*cos(a_z)*(h_z)^2).^2);
c2 = 1 + ((lambda_x)*(lambda_y)*(lambda_z)*(h_x)^2*(h_y)^2*(h_z)^2) - ...
(1/2)*((lambda_x)*cos(a_x)*(h_x)^2).*((lambda_y)*cos(a_y)*(h_y)^2) ...
.*((lambda_z)*cos(a_z)*(h_z)^2) - (1/2)*(((lambda_x)*cos(a_x)*(h_x)^2).^2)...
.*(((lambda_y)*cos(a_y)*(h_y)^2).^2).*(((lambda_z)*cos(a_z)*(h_z)^2).^2);
c3 = 1 + ((lambda_x)*(lambda_y)*(lambda_z)*(h_x)^2*(h_y)^2*(h_z)^2) - ...
(1/2)*((lambda_x)*cos(a_x)*(h_x)^2).*((lambda_y)*cos(a_y)*(h_y)^2) ...
.*((lambda_z)*cos(a_z)*(h_z)^2) - (1/2)*(((lambda_x)*cos(a_x)*(h_x)^2).^2)...
.*(((lambda_y)*cos(a_y)*(h_y)^2).^2).*(((lambda_z)*cos(a_z)*(h_z)^2).^2);
% went out to length 2m
Z1 = sqrt(1./(2*m_y)).*real(fft(V1./c1, 2*m_y,1));
Z1 = Z1([1:2*m_y 1],:,:);
Z2 = sqrt(1./(2*m_y)).*real(fft(V2./c2, 2*m_y,1));
Z2 = Z2([1:2*m_y 1],:,:);
Z3 = sqrt(1./(2*m_y)).*real(fft(V3./c3, 2*m_y,1));
Z3 = Z3([1:2*m_y 1],:,:);
R1 = sqrt(1./(2*m_x)).*real(fft(Z1, 2*m_x,2));
R1 = R1(:,[1:2*m_x 1],:);
R2 = sqrt(1./(2*m_x)).*real(fft(Z2, 2*m_x,2));
R2 = R2(:,[1:2*m_x 1],:);
R3 = sqrt(1./(2*m_x)).*real(fft(Z3, 2*m_x,2));
R3 = R3(:,[1:2*m_x 1],:);
S1 = sqrt(2./m_z).*imag(fft(R1, 2*m_z,3));
S1 = S1(:,:,[1:2*m_z 1]);
S2 = sqrt(2./m_z).*imag(fft(R2, 2*m_z,3));
S2 = S2(:,:,[1:2*m_z 1]);
S3 = sqrt(2./m_z).*imag(fft(R3, 2*m_z,3));
S3 = S3(:,:,[1:2*m_z 1]);
w1(:,:,:,t+1) = S1;
w2(:,:,:,t+1) = S2;
w3(:,:,:,t+1) = S3;
end
%crop w1,w2,w3
w1 = w1(1:m_y+1,1:m_x+1,1:m_z+1,:);
w2 = w2(1:m_y+1,1:m_x+1,1:m_z+1,:);
w3 = w3(1:m_y+1,1:m_x+1,1:m_z+1,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = 10;
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 = [10,10,10];
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);
%compute the 8 corners
y(1,1,1) = (2*x(1,2,1)-2*x(1,1,1))/(h_x.^2)+(2*x(2,1,1)-2*x(1,1,1))/(h_y.^2)...
+(2*x(1,1,2)-2*x(1,1,1))/(h_z.^2);
y(m,1,1) = (2*x(m,2,1)-2*x(m,1,1))/(h_x.^2)+(2*x(m-1,1,1)-2*x(m,1,1))/(h_y.^2)...
+(2*x(m,1,2)-2*x(m,1,1))/(h_z.^2);
y(1,n,1) = (2*x(1,n-1,1)-2*x(1,n,1))/(h_x.^2)+(2*x(2,n,1)-2*x(1,n,1))/(h_y.^2)...
+(2*x(1,n,2)-2*x(1,n,1))/(h_z.^2);
y(1,1,p) = (2*x(1,2,p)-2*x(1,1,p))/(h_x.^2)+(2*x(2,1,p)-2*x(1,1,p))/(h_y.^2)...
+(2*x(1,1,p-1)-2*x(1,1,p))/(h_z.^2);
y(m,n,1) = (2*x(m,n-1,1)-2*x(m,n,1))/(h_x.^2)+(2*x(m-1,n,1)-2*x(m,n,1))/(h_y.^2)...
+(2*x(m,n,2)-2*x(m,n,1))/(h_z.^2);
y(1,n,p) = (2*x(1,n-1,p)-2*x(1,n,p))/(h_x.^2)+(2*x(2,n,p)-2*x(1,n,p))/(h_y.^2)...
+(2*x(1,n,p-1)-2*x(1,n,p))/(h_z.^2);
y(m,1,p) = (2*x(m,2,p)-2*x(m,1,p))/(h_x.^2)+(2*x(m-1,1,p)-2*x(m,1,p))/(h_y.^2)...
+(2*x(m,1,p-1)-2*x(m,1,p))/(h_z.^2);
y(m,n,p) = (2*x(m,n-1,p)-2*x(m,n,p))/(h_x.^2)+(2*x(m-1,n,p)-2*x(m,n,p))/(h_y.^2)...
+(2*x(m,n,p-1)-2*x(m,n,p))/(h_z.^2);
%compute the 12 edges
for i = 2:n-1
y(1,i,1) = (x(1,i-1,1)+x(1,i+1,1)-2*x(1,i,1))/(h_x.^2) + ...
(2*x(2,i,1)-2*x(1,i,1))/(h_y.^2) + (2*x(1,i,2)-2*x(1,i,1))/(h_z.^2);
end
for i = 2:n-1
y(m,i,1) = (x(m,i-1,1)+x(m,i+1,1)-2*x(m,i,1))/(h_x.^2) + ...
(2*x(m-1,i+1,1)-2*x(m,i,1))/(h_y.^2) + (2*x(m,i,2)-2*x(m,i,1))/(h_z.^2);
end
for i = 2:n-1
y(1,i,p) = (x(1,i-1,p)+x(1,i+1,p)-2*x(1,i,p))/(h_x.^2) + ...
(2*x(2,i+1,p)-2*x(1,i,p))/(h_y.^2) + (2*x(1,i,p-1)-2*x(1,i,p))/(h_z.^2);
end
for i = 2:n-1
y(m,i,p) = (x(m,i-1,p)+x(m,i+1,p)-2*x(m,i,p))/(h_x.^2) + ...
(2*x(m-1,i+1,p)-2*x(m,i,p))/(h_y.^2) + (2*x(m,i,p-1)-2*x(m,i,p))/(h_z.^2);
end
for j = 2:m-1
y(j,1,1) = (x(j-1,1,1)+x(j+1,1,1)-2*x(j,1,1))/(h_y.^2) + ...
(2*x(j,2,1)-2*x(j,1,1))/(h_x.^2) + (2*x(j,1,2)-2*x(j,1,1))/(h_z.^2);
end
for j = 2:m-1
y(j,n,1) = (x(j-1,n,1)+x(j+1,n,1)-2*x(j,n,1))/(h_y.^2) + ...
(2*x(j,n-1,1)-2*x(j,n,1))/(h_x.^2) + (2*x(j,n,2)-2*x(j,n,1))/(h_z.^2);
end
for j = 2:m-1
y(j,1,p) = (x(j-1,1,p)+x(j+1,1,p)-2*x(j,1,p))/(h_y.^2) + ...
(2*x(j,2,p)-2*x(j,1,p))/(h_x.^2) + (2*x(j,1,p-1)-2*x(j,1,p))/(h_z.^2);
end
for j = 2:m-1
y(j,n,p) = (x(j-1,n,p)+x(j+1,n,p)-2*x(j,n,p))/(h_y.^2) + ...
(2*x(j,n-1,p)-2*x(j,n,p))/(h_x.^2) + (2*x(j,n,p-1)-2*x(j,n,p))/(h_z.^2);
end
for k = 2:p-1
y(1,1,k) = (x(m,1,k-1)+x(m,1,k+1)-2*x(m,1,k))/(h_z.^2) + ...
(2*x(m,2,k)-2*x(m,1,k))/(h_x.^2) + (2*x(m-1,1,k)-2*x(m,1,k))/(h_y.^2);
end
for k = 2:p-1
y(m,1,k) = (x(1,1,k-1)+x(1,1,k+1)-2*x(1,1,k))/(h_z.^2) + ...
(2*x(1,2,k)-2*x(1,1,k))/(h_x.^2) + (2*x(2,1,k)-2*x(1,1,k))/(h_y.^2);
end
for k = 2:p-1
y(1,n,k) = (x(1,n,k-1)+x(1,n,k+1)-2*x(1,n,k))/(h_z.^2) + ...
(2*x(1,n-1,k)-2*x(1,n,k))/(h_x.^2) + (2*x(2,n,k)-2*x(1,n,k))/(h_y.^2);
end
for k = 2:p-1
y(m,n,k) = (x(m,n,k-1)+x(m,n,k+1)-2*x(m,n,k))/(h_z.^2) + ...
(2*x(m,n-1,k)-2*x(m,n,k))/(h_x.^2) + (2*x(m-1,n,k)-2*x(m,n,k))/(h_y.^2);
end
%compute the inside
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