% infiniteslab.m Last Modified 7/18/2008 % Conduction of heat in an infinite slab with initial constant % temperature T0 on (-L,L) and zero elsewhwere % The diffusivity is kappa=k/(rho*c) % This m-file was written at the University of Wyoming in the Electrical % and Computer Engineering Department and is to be distributed without % cost. clear all set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); xoverL=linspace(0,3,91); sn=[]; nterm=10; ftime=1; ntime=51; tau=linspace(0,ftime,ntime); tau(1,1)=.00001; one=ones(1,91); framedata=[]; for k=1:ntime den=2*sqrt(tau(1,k)); V=0.5*(erf((one-xoverL)/den)+erf((one+xoverL)/den)); framedata=[framedata;V]; end framedata2=[fliplr(framedata(:,2:91)) framedata]; xoverL2=linspace(-3,3,181); figure(1);clf;% Plot T(x,t) vs x for various t xlabel('Dimensionless Distance, x/L') hold on ylabel('Dimensionless Temperature, T(x,t)/T_0') for k=1:2:51 plot(xoverL2,framedata2(k,:),'k') hold on axis([-3 3 -.1 1.2]); hold on end text(-.95,.05,['\kappat/L^2 = ' num2str(0)]) text(-.2 ,.45,['\kappat/L^2 = ' num2str(1)]) set(gca,'Box','on') hold off text(-2.5,1.1,'Press Enter to Continue') pause figure(2);clf;%Plot T(x,t) vs t for various x hold on for j=1:12:91 plot(tau,framedata(:,j)) hold on axis([0 1 0,1]); hold on end text(.4,.76,['x/L = ' num2str(0)]) text(.4,.62,['x/L = ' num2str(0.4)]) text(.4,.52,['x/L = ' num2str(.8)]) text(0.4,0.38,['x/L = ' num2str(1.2)]) xlabel('Dimensionless Time, \kappat/L^2') ylabel('Dimensionless Temperature, T(x,t)/T_0') hold off text(.2,.9,'Press Enter to Continue') pause figure(3);clf%now do animation axis([-3 3,-.1 1.2]) hold on box on xlabel('Dimensionless Distance, x/L') hold on ylabel('Dimensionless Temperature, T(x,t)/T_0') hold on L=plot(xoverL2, framedata2(1,:),'k','EraseMode','xor'); hold on texthandl=text(-2.5, 1.1,'Press Enter to Animate'); pause set(texthandl,'String',' '); for i=2:51 set(L,'Ydata',framedata2(i,:)) pause(.05) end hold on plot(xoverL2,framedata2(1,:),'k') hold on set(texthandl,'String','Press Enter to Continue'); hold off pause figure(5);clf;%3-D Plot of T(x,t) vs x and t [X,Y]=meshgrid(tau,(xoverL2)); mesh(X,Y,framedata2') colormap(cool) text(1,-1,-.4,'Distance, x/L','Rotation',12) text(.1, -3,-.4,'Time, \kappat/L^2','Rotation',-32) zlabel('Dimensionless Temperature, T(x,t)/T_0') view(60,30) axis([0 1 -3 3 -.1 1.2]) text(.1,3,1,'Press Enter to Continue','Rotation',-32) pause