2.1)
Determine the coefficients a to d of the central difference scheme:
where ε is the discretization error. What is the order of accuracy of the scheme? Is the scheme diffusive or dispersive? Assume a uniformly spaced grid ( 
is constant).
is constant).Solution
Expand 
with Taylor series from the point 
with Taylor series from the point 
and substitute into (1)
Next we match the coefficients in front of the derivatives 
on the left and right hand side of equation (3) to obtain an algebraic system for the constants 
. We can write it in matrix-vector form as
on the left and right hand side of equation (3) to obtain an algebraic system for the constants . We can write it in matrix-vector form as
Solving the system in MatLab we obtain
syms dx
A = [ 1 1 1 1
-2 -1 1 2
4 1 1 4
-8 -1 1 8 ] ;
f = [ 0 1/dx 0 0 ]';
abcd = A\f
Substituting into (1) we obtain the finite-difference approximation
into (1) we obtain the finite-difference approximation
In order to express the leading term of the approximation error, ε, we also substitute into (3) and obtain
into (3) and obtain
Thus the scheme is fourth-order accurate and the leading order dicretization error has a dispersive nature (scales with an odd derivative of the solution).
2.2)
Determine the coefficients a and b of the following central finite difference scheme for a NON-UNIFORM grid (
const.):
const.):
What is the order of accuracy of this scheme? Would the order of accuracy change for a uniform spacing? Can you conclude any recommendations for designing a non-uniform finite-difference grids?
Notation: 
Solution
The linear system for a and b then reads
syms dx_im1
A = [ 1 1
-dx_im1 dx ] ;
f = [ 0 1 ]';
ab = A\f
Substituting these coefficients we obtain a centered finite-difference approximation for non-uniform grids
Substituting the Taylor expansion of 
and 
we can express the discretization error
and we can express the discretization error
One can see already that the leading term of the discretization error vanishes if 
. We would therefore conclude that the grid spacing should not increase nor decrease too rapidly in space, otherwise the discretization scheme becomes too diffusive (due to the even derivative in the leading order error term). To see the order of accuracy more clearly we define a growth factor 
to obtain
. We would therefore conclude that the grid spacing should not increase nor decrease too rapidly in space, otherwise the discretization scheme becomes too diffusive (due to the even derivative in the leading order error term). To see the order of accuracy more clearly we define a growth factor to obtain
The scheme is therefore first order accurate for non-uniform grids and second order accurate for uniform grids.
2.3)
For the function 
compute the first derivative 
at 
using the schemes below with 
:
compute the first derivative at using the schemes below with :
and compare the results with the exact solution. For each scheme, plot the dependence of the discretization error on in logarithmic scale. Can you determine the order of accuracy of the scheme from the plot?
in logarithmic scale. Can you determine the order of accuracy of the scheme from the plot?Solution
u = @(x) sin(pi*x);
dudx = @(x) pi*cos(pi*x);
x = 0.4;
minpow = -5;
maxpow = -1;
resol = 20;
dx = logspace( minpow, maxpow, resol );
dudx1= ( u(x+dx)- u(x ) ) ./ dx ;
dudx2= ( u(x+dx)- u(x-dx) ) ./ ( 2*dx) ;
dudx3= (u(x-2*dx)-8*u(x-dx)+8*u(x+dx)-u(x+2*dx)) ./ (12*dx) ;
e1 = abs( dudx1 - dudx(x) ) ;
e2 = abs( dudx2 - dudx(x) ) ;
e3 = abs( dudx3 - dudx(x) ) ;
clf;
loglog( dx,e1,'.', dx,e2,'.', dx,e3,'.', dx,dx.^1,'k-', dx,dx.^2,'k-', dx,dx.^4,'k-' );
xlabel('$\Delta x$', 'Interpreter','latex');
ylabel('$\epsilon$', 'Interpreter','latex');
legend('(a)','(b)','(c)', 'Location','southeast');
annotation('textbox',[.3 .06 .2 .2],'String','$\Delta x^4$', 'Interpreter','latex','FitBoxToText','on','EdgeColor','none');
annotation('textbox',[.3 .38 .2 .2],'String','$\Delta x^2$', 'Interpreter','latex','FitBoxToText','on','EdgeColor','none');
annotation('textbox',[.3 .55 .2 .2],'String','$\Delta x$ ', 'Interpreter','latex','FitBoxToText','on','EdgeColor','none');
Note that the slope of error ε with reveals the order of accuracy of the scheme.
reveals the order of accuracy of the scheme.2.4) Steady heat conduction
The steady heat conduction in one dimension is governed by:
where T is temperature, α is thermal diffusivity and q is internal heat generation. With finite differences we can discretize the equation as
Let us assume the case 
with fixed temperature on boundaries
with fixed temperature on boundaries
We write the discrete problem in matrix-vector form
and solve it with MatLab. First we need an analytical solution
syms Ta(x)
F = sin(2*pi*x); % Source term
x1 = 0; xN = 1; % Boundaries of the domain
Ta = dsolve( -diff(Ta,2)==F, [Ta(x1)==0,Ta(xN)==0] ) % Analytical solution
figure(1); clf; hold on; % Clean the figure window
fplot( Ta, [x1,xN], 'DisplayName', 'analytical' ); % Plot of the analytical solution
Ta = matlabFunction(Ta); % Convert symbolics to function handles
F = matlabFunction(F); % otherwise the evaluation takes too long
in order to compute the discretization error
N = [5,10,20,40]; % Number of grid points
eps= NaN(size(N)); % Numerical error
for k=1:length(N)
dx = (xN-x1) / (N(k)-1); % Grid spacing
x = (x1:dx:xN)'; % Grid points
e = ones(N(k),1); % Auxiliary column vector of N ones
A = spdiags([e -2*e e]/dx^2,[-1 0 1],N(k),N(k)); % Matrix of discrete second derivative
f = F(x); % Vector of right hand sides
A(1 ,:) = 0; A(1 ,1 ) = 1; f(1 ) = 0; % homogeneous Dirichlet BC at x1
A(end,:) = 0; A(end,end) = 1; f(end) = 0; % homogeneous Dirichlet BC at xN
T = -A\f; % solution T = A^-1 * f
plot(x, T,'.','DisplayName',num2str(N(k),'%1.1e')); % plot discrete solution
eps(k) = max( abs(T-Ta(x)) ) / max(Ta(x)); % discretization error
end
xlabel('x'); ylabel('T'); legend('show');
figure(2); clf;
loglog(N,eps,'-*', N,0.01*ones(size(N)),'k--', N,N.^(-2),'k-');
xlabel('$N$','Interpreter','latex'); ylabel('$\epsilon$','Interpreter','latex');
grid on;
legend({'Discretization error','$\epsilon=1\%$','$N^{-2}$'},'Interpreter','latex')
From the second figure it is clear that we need at least 20 grid points to reduce the discretization error below 1%. We can also see from the slope of error decay that the scheme is 2nd order accurate.
Let us now replace the source term
syms Ta(x)
F = @(x) ones(size(x)); % Source term
Ta = dsolve( -diff(Ta,2)==F, [Ta(x1)==0,Ta(xN)==0] ) % Analytical solution
figure(1); clf; hold on; % Clean the figure window
fplot( Ta, [x1,xN], 'DisplayName', 'analytical' ); % Plot of the analytical solution
Ta = matlabFunction(Ta); % Convert symbolics to function handles
such that the analytical solution is a parabola. Then compute the discretization error as before
N = [3,6,12]; % Number of grid points
eps= NaN(size(N)); % Numerical error
for k=1:length(N)
dx = (xN-x1) / (N(k)-1); % Grid spacing
x = (x1:dx:xN)'; % Grid points
e = ones(N(k),1); % Auxiliary column vector of N ones
A = spdiags([e -2*e e]/dx^2,[-1 0 1],N(k),N(k)); % Matrix of discrete second derivative
f = F(x); % Vector of right hand sides
A(1 ,:) = 0; A(1 ,1 ) = 1; f(1 ) = 0; % homogeneous Dirichlet BC at x1
A(end,:) = 0; A(end,end) = 1; f(end) = 0; % homogeneous Dirichlet BC at xN
T = -A\f; % solution T = A^-1 * f
plot(x, T,'.','DisplayName',num2str(N(k),'%1.1e')); % plot discrete solution
eps(k) = max( abs(T-Ta(x)) ) / max(Ta(x)); % discretization error
end
xlabel('x'); ylabel('T'); legend('show');
figure(2); clf;
loglog(N,eps,'-*');
xlabel('$N$','Interpreter','latex'); ylabel('$\epsilon$','Interpreter','latex');
Note that no matter how many grid points we take, the numerical solution always coincides with the analytical solution. Even with only 1 interior node, the discretization error drops below the double precision. One can show using Taylor expansion that the error of this discretization scheme scales with 3rd and higher derivatives of the solution, which are zero in this case.
Important Note: Interpolating the numerical solution (e.g. for plotting purposes) linearly on a small number of grid points can create additional errors, although the numerical solution is exact in this case. Always interpolate numerical solutions on coarse grids with an interpolant of the same order as the discretization scheme!
2.5) Friedrich's equation
a) Assume that f is a passive scalar. Explain the physical meaning of each term.
Friedrich's equation is equivalent to a classical convection-diffusion problem
- diffusion:
- convection: -
- generation: -a
b) Discretize the equation using finite differences on a uniform grid
Common choice for finite differences are second-order centered schemes. Matrix representation of the discrete second derivative was given in the exercise 2.4) as
Second order discretization of first derivative was given in 2.3b) as
epsilon = 0.0001;
N = 101;
dx = 1/(N-1);
x = (0:dx:1)';
e = ones(N,1);
a = 1/dx^2 * e; % lower diagonal
b = -2/dx^2 * e; % main diagonal
c = 1/dx^2 * e; % upper diagonal
D2 = spdiags( [a b c], [-1 0 1], N, N);
m = -1/(2*dx)*e;
n = 1/(2*dx)*e;
D1 = spdiags( [m n], [-1 1], N, N);
A = epsilon*D2 + D1;
f = 2 .* e;
A(1,1) = 1;
A(1,2) = 0;
f(1) = 0;
A(end,end ) = 1;
A(end,end-1) = 0;
f(end) = 1;
T = A\f;
clf; plot(x,T);
