clear all; %close all; % Set non-dimensional thermal coefficient k = 1.0; % this is really k/(rho*cp) % Set length of bar L = 1.0; % non-dimensional % Set initial temperature Tinit = 400; % Set left and right temperatures for t>0 Tleft = 800; Tright = 1000; % Set up grid size Nx = 1001; h = L/Nx; x = linspace(0,L,Nx+1); % Set timestep size dt = 0.001; % Calculate number of iterations (Nmax) needed to iterate to t=Tmax Tmax = 0.5; Nmax = floor(Tmax/dt); t = linspace(0,Nmax*dt,Nmax+1); Tmid = zeros(Nmax,1); Tmid(1) = Tinit; % Initialize a sparse matrix to hold stiffness & identity matrix A = spalloc(Nx-1,Nx-1,3*(Nx-1)); I = speye(Nx-1); % These commands allocate full (non-sparse) matrices %A = zeros(Nx-1,Nx-1); %I = eye(Nx-1); % Calculate stiffness matrix for ii = 1:Nx-1, if (ii > 1), A(ii,ii-1) = k/h^2; end if (ii < Nx-1), A(ii,ii+1) = k/h^2; end A(ii,ii) = -2*k/h^2; end % Set forcing vector b = zeros(Nx-1,1); b(1) = k*Tleft/h^2; b(Nx-1) = k*Tright/h^2; % Set initial vector v = Tinit*ones(Nx-1,1); figure(1); % Iterate in time for n = 1:Nmax, % Backward Euler % v = (I - dt*A)\(v + dt*b); % Forward Euler %v = v + dt*(A*v+b); % Trapezoidal rhs = v + 0.5*dt*A*v + dt*b; G = I - 0.5*dt*A; v = G\rhs; % 2nd order bd % if (n==1), % v1 = v; % v = (I - dt*A)\(v1 + dt*b); % else, % v0 = v1; % v1 = v; % v = (I - 2/3*dt*A)\(4/3*v1 - 1/3*v0 + 2/3*dt*b); % end T = [Tleft; v; Tright]; plot(x,T,'*'); axis([0,1,400,1200]); drawnow; Tmid(n+1) = T(floor(0.5*(Nx+1))); end figure(2); plot(t,Tmid);hold on; figure(3); plot(x,T); hold on; axis([0,1,400,1200]);