% Matlab program heat4 % % This program integrates the heat equation u_t - ku_xx = q(x) on % the interval [0, 10] with Neumann boundary conditions u_x(0,t) = % u_x(10,t) = 0. The Crank-Nicholson method is used with delx = .05 % and delt = .05. % At run time user must enter the value of the diffusion % coefficient k. User also must enter the number of time steps % n1, n2, n3, and n4 between snapshots as in program heat3. % Output is written into the vectors 'snap0',snap1', 'snap2','snap3' % and 'snap4'. The vectors may be plotted alone or together. % The initial data comes from the file f.m. User must % provide a function mfile q.m for the source. disp(' ') k = input('enter the value of the diffusion coefficient k ') delt = .05 disp(' Enter the number of time steps between snapshots, n1,n2, n3, n4 ') n = input(' in the form [n1 n2 n3 n4] ') n1 = n(1); n2 = n(2); n3 = n(3); n4 = n(4); t1 = n1*delt t2 = (n1 +n2)*delt t3 = (n1 + n2 +n3)*delt t4 = (n1 + n2 + n3 +n4)*delt delx = .05; r = .5*k*delt/(delx^2); J = 200; x = 0:delx:10; snap0= f(x); qq = q(x); v= snap0; D= sparse(1:J+1, 1:J+1, (1+2*r), J+1,J+1, J+1); D(1,1) = D(1,1) - r; D(J+1, J+1) = D(J+1, J+1) - r; E = sparse(2:J+1, 1:J, -r, J+1,J+1,J+1); A = D + E+ E'; [L,U] = lu(A); for n = 1:n1 b = 2*v + delt*qq; y = L\b'; z = (U\y)'; v = z - v; end snap1 = v; disp('Computed up to time t1 ') for n = n1+1: n1+n2 b = 2*v +delt*qq; y = L\b'; z = (U\y)'; v = z-v; end snap2 = v; disp('Computed up to time t2') for n = n1 + n2+1:n1 + n2 +n3 b = 2*v +delt*qq; y = L\b'; z = (U\y)'; v = z-v; end snap3 = v; disp('Computed up to time t3') for n = n1 + n2 + n3+1:n1 + n2 + n3 +n4 b = 2*v +delt*qq; y = L\b'; z = (U\y)'; v = z-v; end snap4 = v; disp('Computed up to time t4') plot(x,snap0,x,snap1,x,snap2,x,snap3,x,snap4) legend('t0','t1','t2','t3','t4')