2 Consider the function on the interval 22 with h 025 Use th

2. Consider the function: on the interval (-2,2 with h 0.25. Use the forward, backward, and centered finite difference approx- imations for the first and second derivatives so as to graphically illustrate which approximation is most accurate. Graph all three first-derivative finite difference approximations along with the theoretical, and do the same for the second derivative as well.

Solution

clear, clf
C = 1; D = .01; K = .1; A=[]; n = 4;

for i = 1:8
flops(0)
h = 1/n;
%Construct the Tridiagonal System
l=ones(1,n-2); d=-(2+K*h^2/D)*ones(1,n-1); b=[zeros(1,n-2), -C];
%Solve the Tridiagonal System
U = tridiag(l,d,l,b);
p=flops;
%Plot the Approx. Solution
X = 0:h:1; U = [0, U, C];
plot(X,U,\'b\')
grid on, xlabel(\'x-axis\'), ylabel(\'u-axis\'), title(\'Example 2.6\')
hold on
%Save values at the midpoint for error analysis
A=[A;n,p,X(n/2+1),U(n/2+1)];
n = n*2;
end

hold off
%error analysis
[nn,mm]=size(A);
diff(1)=0;ratio(1)=0; ratio(2)=0;

for i=1:nn-1,
diff(i+1)=A(i+1,4)-A(i,4);
end

for i=2:nn-1,
ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,4) + diff(i)/(1-ratio(i+1));
end

z=[A(:,1)\';A(:,2)\';A(:,4)\'; diff; ratio; xacc];
fprintf(1,\'\ \ n flops y diff ratio y acc\ \ \')
fprintf(1,\'%3.0f %7.0f %16.14f %13.10f %13.10f %16.14f\ \',z)

clear, clf
C = 1; D = .01; K = .1; A=[];n = 4;
for i = 1:8
h = 1/n;
flops(0)
%Using Backward Diff at endpoint
%l=[ones(1,n-2),-1]; u=ones(1,n-1); d=[-(2+K*h^2/D)*ones(1,n-1),1];
%b=[zeros(1,n-1), h*C/D];Y = tridiag(l,d,u,b);
%Using an imaginary node at x=1+h
l=[ones(1,n-2),2]; u=ones(1,n-1); d=[-(2+K*h^2/D)*ones(1,n)];
b=[zeros(1,n-1), -2*h*C/D];Y = tridiag(l,d,u,b);
%Using a second order approx at x=1
%l=[ones(1,n-2),-4];u=ones(1,n-1);d=[-(2+K*h^2/D)*ones(1,n-1),3];
%b=[zeros(1,n-1), 2*h*C/D];
%M=zeros(n,n)+diag(d)+diag(l,-1)+diag(u,1);M(n,n-2)=1;
%Y=(M\\b\')\';
p=flops;
X = 0:h:1; Y = [0, Y];
plot(X,Y,\'b\')
grid on, xlabel(\'x-axis\'), ylabel(\'u-axis\'), title(\'Example 2.7\')
hold on
%Save values at the midpoint for error analysis
A=[A;n,p,X(n/2+1),Y(n/2+1)];
n = n*2;
end
hold off
A;
[nn,mm]=size(A);
diff(1)=0;ratio(1)=0; ratio(2)=0;
for i=1:nn-1,
diff(i+1)=A(i+1,4)-A(i,4);
end
for i=2:nn-1,
ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,4) + diff(i)/(1-ratio(i+1));
end
z=[A(:,1)\';A(:,2)\';A(:,4)\'; diff; ratio; xacc];
fprintf(1,\'\ \ n flops u(0.5) diff ratio u acc\ \ \')
fprintf(1,\'%3.0f %7.0f %16.14f %13.10f %13.10f %16.14f\ \',z)