% expm10.m - for Table 7 % This computes sig_N(t) by discretizing the interval (0,a) close('all') clear all % Set up basic data %h=number of divisions per unit length h=20 %a=length of interval a=10 %b=diffusion paramter b=20 % number of points altogether N = h*a-1 M= 2*h*a-1 % set up tridiagonal matrix for N and compute function u=ones(N,1); v=ones(N-1,1); A=-2*(h^2)*b^(-1)*diag(u)+(h^2)*b^(-1)*diag(v,1)+(h^2)*b^(-1)*diag(v,-1); B=(-h/2)*diag(v,1)+(h/2)*diag(v,-1); C=A+B; ft=zeros(N,10); for t=1:10 ft(:,t)=expm((t/10)*C)*ones(N,1); end; % set up tridiagonal matrix for M and compute function h2=2*h; u2=ones(M,1); v2=ones(M-1,1); A2=-2*(h2^2)*b^(-1)*diag(u2)+(h2^2)*b^(-1)*diag(v2,1)+(h2^2)*b^(-1)*diag(v2,-1); B2=(-h2/2)*diag(v2,1)+(h2/2)*diag(v2,-1); C2=A2+B2; gt=zeros(M,10); for t=1:10 gt(:,t)=expm((t/10)*C2)*ones(M,1); end; % compute difference of two functions ht=zeros(N,10); for r=1:N for t=1:10 ht(r,t)=ft(r,t)-gt(2*r,t); end; end; % compute norm of difference sigN=zeros(1,10); for t=1:10 sigN(1,t)=norm(ht(:,t),inf); end; format long sigN' hold on for t=1:10 plot(h^(-1)*(1:N),abs(ht(:,t))); end;