% This code solve the one-dimensional heat diffusion equation % for the problem of a bar which is initially at T=Tinit and % suddenly the temperatures at the left and right change to % Tleft and Tright. % % Upon discretization in space by a finite difference method, % the result is a system of ODE's of the form, % % u_t = Au + b % % The code calculates A and b. Then, these are used in % various different time integration methods to evolve the % temperature in time. See lecture notes for more info. % % Using this code, we will talk about the following issues: % % * Use of explicit methods for stiff systems % * Implementation of implicit methods for linear systems % * Sparse vs. full matrices % * Behavior of trapezoidal when lambda dt -> -infinity. 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 = 1000; h = L/Nx; x = linspace(0,L,Nx+1); % Set timestep size dt = 1e-1; % Calculate number of iterations (Nmax) needed to iterate to t=Tmax Tmax = 0.5; Nmax = floor(Tmax/dt); % 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); % 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; T = [Tleft; v; Tright]; plot(x,T,'*'); axis([0,1,400,1200]); drawnow; end