Skip to content

Commit

Permalink
1
Browse files Browse the repository at this point in the history
  • Loading branch information
Leahy committed Dec 12, 2017
1 parent 8354c38 commit efc768c
Showing 1 changed file with 25 additions and 88 deletions.
113 changes: 25 additions & 88 deletions SE_diff.m
Original file line number Diff line number Diff line change
@@ -1,91 +1,28 @@
function [pw_se,w]=SE_diff(T,P,n)
w = membrane_solution31(T,P,n);
A = zeros(n,n);
h = zeros(n,1);
m = n+1; %Represents the number of nodes

%Top Left Corner Interior Node
for i = 2+m
A(i,i)=1;
A(i,i-1)=0;
A(i,i+m)=1;
A(i,i+m-1)=0;
h(i) = 10;
end
%Top Right Corner Interior Node
for i = m*2-1
A(i,i)=1;
A(i,i-1)=1;
A(i,i+m)=1;
A(i,i+m-1)=1;
h(i) = 10;
end
%Bottom Left Corner Interior Node
for i = m*m-(2*m)+1
A(i,i)=1;
A(i,i-1)=0;
A(i,i+m)=0;
A(i,i+m-1)=0;
h(i) = 10;
end
%Bottom Right Corner Interior Node
for i = (m*m)-m-1
A(i,i)=1;
A(i,i-1)=1;
A(i,i+m)=0;
A(i,i+m-1)=0;
h(i) = 10;
end

%The following represent the values for the interior nodes

%Top Interior Row
for i = 3+m:m*2-2
A(i,i)=1;
A(i,i-1)=1;
A(i,i+m)=1;
A(i,i+m-1)=1;
h(i) = 10;
end
%Bottom Interior Row
for i = m*m-(2*m)+2:(m*m)-m-2
A(i,i)=1;
A(i,i-1)=1;
A(i,i+m)=0;
A(i,i+m-1)=0;
h(i) = 10;
end
%Left Interior Column
for i = (2*m)+2,(2*m)+2:(m*m)-(2*m)+2
A(i,i)=1;
A(i,i-1)=0;
A(i,i+m)=1;
A(i,i+m-1)=0;
h(i) = 10;
end
%Right Interior Column
for i = (3*m)-1,(2*m)+2:(m*m)-(2*m)+2
A(i,i)=1;
A(i,i-1)=1;
A(i,i+m)=1;
A(i,i+m-1)=1;
h(i) = 10;
end
%Inerior Nodes
for i = (2*m)+3:(3*m)-2,(2*m)+3:(m*m)-(3*m)+3
A(i,i)=1;
A(i,i-1)=1;
A(i,i+m)=1;
A(i,i+m-1)=1;
h(i) = 10;
end


dw_dx = (A(i,i)-A(i-1)+A(i,i+m)-A(i,i+m-1))/(2.*h);
dw_dx'

end
function [pw_se,w] = SE_diff(T,P,n)

% Call function from part c to obtain vector w
w = membrane_solution(T,P,n);
% Create All Nodes Matrix (Interior & Exterior)
m = zeros(n+2);
c = vec2mat(w,n); % Vector w rearranged to n x n matrix.
m(2:n+1,2:n+1) = c; % insert matrix c into middle of m matrix
disp(m); % symmetric matrix

h = 10/(n+1);
% Compute dw_dx
dx = zeros(n+2,n+1);
for i = 2:n+2
dx(:,i-1)= (m(:,i)-m(:,i-1))/h;
end
disp(dx)
dw_dx = zeros(n+1,n+1);
for i = 1:n+1
%dw_dx((i:n+1),i) = (dx(i,i)+dx(i+1,i))/2;
%dw_dx((i-1:n+1),i+1) = (dx(i,2)+dx(i+1,2))/2;
%dw_dx(i,3) = (dx(i,3)+dx(i+1,3))/2;
%dw_dx(i,4) = (dx(i,4)+dx(i+1,4))/2;
%i = i+1;
end
disp(dw_dx)


% Compute dw_dy

0 comments on commit efc768c

Please sign in to comment.