Lecture1. Things to do in class.

Contents

Number 1

clear all; close all; % clear all to start fresh
% Note: n equal spacing means n+1 points
% as an example let use use a = -1 and b = 1
a = -1; b = 1;
n = 10; h = (b-a)/n;
one = ones(n+1,1); % create a column vector of size n+1 by 1 of ones
D = spdiags([-one one],[0 1],n+1,n+1);
D(end,end-1) = -1; D(end,end) = 1; % Modify last row (one-sided weights)
D = (1/h)*D; % Don't forget the 1/h factor
spy(D) % Plot the structure of the sparse matrix

Number 2

format long e % format number in long format
u = @(x) exp(-x.^2); % Define the function
ux = @(x) -2*x.*exp(-x.^2); % Define the exact derivative of the function
x = linspace(a,b,n+1).'; % Create n+1 equally-spaced points in [a,b]
uux = D*u(x); % Compute numerical derivative
err = norm(uux-ux(x),Inf) % Compute infnorm (u'_num - u'_exact)
figure
plot(x,uux,'-*r',x,ux(x),'-b')

err =

     1.960528042383842e-01

Number 3

Can you check what is the slope of the error convergence ?

n = 10.^(1:6); % Create n = 10,100,1000,10000,...
h = (b-a)./n;
err = zeros(length(n),1); % Pre-allocate memory to store errors
for i=1:length(n)
    one = ones(n(i)+1,1);
    x = linspace(a,b,n(i)+1).';
    D = spdiags([-one one],[0 1],n(i)+1,n(i)+1);
    D(end,end-1) = -1; D(end,end) = 1; D = (1/h(i))*D;
    uux = D*u(x);
    err(i) = norm(uux-ux(x),Inf);
end
figure
loglog(h,err,'-o') % Plot the error in logscale
xlabel('h'); ylabel('Error')

Number 4

In periodic case, x0 and xn are the same point. So, you can remove either x0 or xn

u = @(x) sin(pi*x);
ux = @(x) pi*cos(pi*x);
n = 100; h = (b-a)/n;
one = ones(n,1); % create a column vector of size n+1 by 1 of ones
D = spdiags([-one one],[0 1],n,n);
D(end,1) = 1; % Wrap around the matrix
D = (1/h)*D; % Don't forget the 1/h factor
figure
spy(D) % Plot the structure of the sparse matrix
x = linspace(a,b,n+1).'; % Create n+1 equally-spaced points in [a,b]
uux = D*u(x(1:n)); % Compute numerical derivative
err = norm(uux-ux(x(1:n)),Inf) % Compute infnorm (u'_num - u'_exact)
figure
plot(x(1:n),uux,'-*r',x(1:n),ux(x(1:n)),'-b')

err =

     9.866357858642535e-02

Check the error convergence of number 4

n = 10.^(1:6); % Create n = 10,100,1000,10000,...
h = (b-a)./n;
err = zeros(length(n),1); % Pre-allocate memory to store errors
for i=1:length(n)
    one = ones(n(i),1);
    x = linspace(a,b,n(i)+1).';
    D = spdiags([-one one],[0 1],n(i),n(i));
    D(end,1) = 1; D = (1/h(i))*D;
    uux = D*u(x(1:n(i)));
    err(i) = norm(uux-ux(x(1:n(i))),Inf);
end
figure
loglog(h,err,'-o') % Plot the error in logscale
xlabel('h'); ylabel('Error')

Published with MATLAB® R2014b