% program heat1 % % This program plots a Riemann sum which approximates the % convolution of the heat kernal with a function ff(x). The % Riemann sum uses the midpoint rule. % The initial data is written into the vector snap0 and the % approximate solution at time t is written into the vector snap1. % At run time, user must enter the diffusion coefficient k, % the time t > 0 of the snapshot, and the number of terms n (must be even) % in the Riemann sum. The program plots snap0 together with % snap1 on the interval [-10, 10]. % The program also computes the Gaussian approximation at time t. % To compare the solution snap1 with the Gaussian approximation % use plot(x,snap1,x,gauss). % % There are two choices of initial data. The first choice is % a step function % % f(x) = -1 for -1 < x < 0 % f(x) = 2 for 0 < x < 1 % This function is "built in" to the program. If you wish to % have other initial data, you must provide an mfile ff.m or an % inline function ff = inline . . . . % Even if the function ff is defined for all x, the % program heat1 only uses the values of ff for -1 < x < 1, so that it is % equivalent to setting the function to zero outside the interval % [-1,1]. disp( ' ') k = input('enter the value of the diffusion constant k ') t = input('enter the value of the time t to view solution ') n = input('enter the number of points y(j) ') dely = 2.0/n; delx = .01; x = -10:delx:10; y = -1 +.5*dely :dely :1 -.5*dely; disp(' Is the initial data the given step function ? ') mm = input('Enter 1 if yes, 0 if not '); if mm == 1 snap0 = (x < -1) -(x < 0) + 2*( ( x < 1) - (x < 0) ); s = zeros(size(x)); for i = 1:(n/2) s = s - exp(-(x-y(i)).^2/(4*k*t)); end for i = n/2 +1 : n s = s + 2*exp(-(x-y(i)).^2/(4*k*t)); end snap1 = (1.0/sqrt(4*pi*k*t))*s*dely; Q = 1 mean = 1.5 p0= - 1.9167 gauss = Q*(1-p0/(4*k*t))*exp(-(x-mean).^2/(4*k*t) )/sqrt(4*pi*k*t); else snap0 = ( (x < 1) - (x < -1) ).*ff(x); s = zeros(size(x)); for i = 1:n s = s + exp(-(x-y(i)).^2/(4*k*t))*ff(y(i)); end snap1 = (1.0/sqrt(4*k*pi*t))*s*dely; simp = 2*ones(1,201); simp(2:2:200) = 4*ones(100,1); simp(1) = 1; simp(201) = 1; xx = linspace(-1,1, 201)'; funval = ff(xx); Q = (1/300)*simp*funval if abs(Q) < eps disp(' The quantity Q is too close to zero; program cannot calculate the gaussian approximation ') plot(x,snap0, x,snap1) else mean = (1/300)*simp*(xx.*funval)/Q p0 = (1/300)*simp*(xx.^2.*funval) - mean^2 gauss = Q*(1-p0/(4*k*t))*exp(-(x-mean).^2/(4*k*t) )/sqrt(4*pi*k*t); end end plot(x,snap0,x,snap1)