findtemp.m

solvin Parbolic forms of PDE explicitly. The error is O( h2+k2) or O(2,2) , it doesn鈥檛 develop the

function [T,x]=findtemp(a)
%a=3.5;

% input variables
%TL=Left hadn boundary condition
%TR=right hand boundary condition
%TA=ambient temperature
%n=number of internal nodes
%ns=number of segments

 TL = 1; TR=0 ; TA =20; n = 100;
 ns = n + 1;L=1;
 %grid
 dx=L/ns;  x = [1:n]*dx;
 % create A and B matrix
B = zeros(n,1); B(1) = TL ; B(n) = -TR ;           %what is '-' for?
A = zeros(n); A(1,1) = 2+a*dx^2; A(1,2) = -1;
for i=2:n-1
A(i,i-1) = -1;
A(i,i) = 2+a*dx^2;
A(i,i+1) = -1;
end
A(n,n-1) = -1; A(n,n) = 2+a*dx^2;
T= gauss(A,B);
%disp(T); disp(x);
plot(x,T,'ok')
xlabel('distance x')
ylabel('Temperature T ')
legend('Numerical')
end


function x = gauss(A,B)
% gauss: gauss-elimination with partial pivoting

% x = gauss(A,B): gauss elimination with partial piuvoting
% input:
% A = coefficient matrix
% B = right hand side vector
% intermediate:
% m = number of row/equations
% n = number of column/variables in the coefficient matrix
% nb = dimension of the augmented matrix
% Aug = augmented matrix
% output:
% x = solution vector
%
[m,n] = size(A);
if m~=n, error('Matrix A must be square'); end 
nb = n+1;
Aug = [A B];
% forward elimination
for k = 1:n-1
%partial pivoting
[big,i] = max(abs(Aug(k:n,k)));
ipr = i + k - 1;
if ipr~=k
Aug([k,ipr],:)=Aug([ipr,k],:);
end
for i=k+1:n
factor = Aug(i,k)/Aug(k,k);
Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);
end
end
% back substitution
x = zeros(n,1); x(n) = Aug(n,nb)/Aug(n,n);
for i=n-1:-1:1
x(i) =(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
end
end