% Matlab program heat3 % % This program integrates the heat equation u_t -ku_xx = 0 % on the interval [0,10] using the Crank-Nicholson method. % The implicit set of equations are solved at each time step % using an LU factorization. The time step delt = .05, and % delx = .05 as well. % The program makes a snapshot, 'snap0', of the initial data % and four later snapshots 'snap1','snap1','snap3,'snap4' at % times t1,t2,t3, and t4. You may plot the vectors snap0, snap1, etc., % separately or together. % At run time the user must enter the number n1 of time steps % to the time t1 of the first snapshot, snap1, the number n2 of % time steps from time t1 to time t2 of the second snapshot, etc., up to % n4. Example: If n1 = 20, n2 = 20, n3 = 10, n4 = 15, the snap- % shots are taken at times t1 = 1, t2 = 2, t3 = 2.5, and t4 = 3.25. % User must provide a function file f.m for the initial data % f(x). The function must be array smart. % User must also provide function files left.m and right.m for the % values u(0,t) and u(10,t). The boundary data functions do % not need to be array smart. disp(' ') k = input( 'enter the value of the diffusion coefficient k ') delt = .05 disp(' Enter the number of timesteps between snapshots ') 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 = t1 +n2*delt t3 = t2 +n3*delt t4 = t3 +n4*delt delx = .05 r = .5*k*delt/(delx^2); J= 10/delx; x = 0:delx:10; snap0 = f(x); v = snap0(2:1:J); D = sparse(1:J-1, 1:J-1,(1+2*r), J-1, J-1, J-1); E = sparse(2:J-1, 1:J-2, -r, J-1, J-1, J-1); A = D + E +E'; [L,U] = lu(A); % Note: the original index j runs from j = 1 ( x = 0) to j = J ( x = 10 -delx). % The index in v (and hence in b and z) runs from j = 1 (x= delx) to % j = J-1 (x = 10-delx). for n = 1:n1 b = 2*v ; b(1) = b(1) + r*( left(n*delt) +left((n-1)*delt) ); b(J-1) = b(J-1) + r*( right(n*delt) +right((n-1)*delt) ) ; y = L\b'; z = (U\y)'; v = z- v; end snap1 = [left(n*delt),v,right(n*delt)]; disp('Computed up to time t1') for n = n1+1: n1+n2 b = 2*v; b(1) = b(1) + r*( left(n*delt) + left((n+1)*delt) ); b(J-1) = b(J-1) + r*( right(n*delt) + right((n+1)*delt) ); y = L\b'; z = (U\y)'; v = z-v; end disp('Computed up to time t2 ') snap2 = [left(n*delt),v,right(n*delt)]; for n = n1+n2+1: n1 + n2 +n3 b = 2*v; b(1) = b(1) +r*( left(n*delt) + left((n+1)*delt) ); b(J-1) = b(J-1) + r*( right(n*delt) + right((n+1)*delt) ); y = L\b'; z = (U\y)'; v = z-v; end disp('Computed up to time t3 ') snap3 = [left(n*delt),v,right(n*delt)]; for n = n1 + n2 +n3+1: n1+n2+n3+n4 b = 2*v; b(1) = b(1) +r*( left(n*delt) + left((n+1)*delt) ); b(J-1) = b(J-1) + r*( right(n*delt) + right((n+1)*delt) ); y = L\b'; z = (U\y)'; v = z-v; end disp('Computed up to time t4 ') snap4 = [left(n*delt),v,right(n*delt)]; plot(x,snap0,x,snap1,x,snap2,x,snap3,x,snap4)