%
% FD.M Example of Finite Difference Method for Solving ODEs
%
% This demo uses the finite difference technique as the solution method
% for a particular ODE. The mesh interval is constant. Two different
% approximations for the first derivative are used. This comparison
% demonstrates that the best approximation for the derivative terms
% leads to a more efficient numerical scheme. One should always use
% the best approximation available for the given problem!
%
% The sensitivity to mesh size is also addressed by looping over four
% different values of N (number of unknowns in the problem).
%
% Description of given differential equation
% y'' + 3xy' + 7y = cos(2x) with y(0) = 1 and y(pi) = 0
%
% A related Matlab file, FDS.M, is also available (this contains a simple
% version of this example - Case A only - with no mesh sensitivity study).
%
% File prepared by J. R. White, UMass-Lowell (Aug. 2003).
%
%
% getting started
clear all, close all, nfig = 0;
%
% set domain limits and boundary conditions
xo = 0; xf = pi; yxo = 1; yxf = 0;
%
% loop over four values of N to evaluate sensitivity to mesh size
NN = [20 40 80 160];
%
for n = 1:4
N = NN(n);
%
% compute interval size and discrete x vector
dx = (xf-xo)/(N+1); dx2 = dx*dx; x = (xo+dx):dx:(xf-dx);
%
%
% Approximations: y'' = [y(i-1) - 2y(i) + y(i+1)]/(dx^2) error => (dx^2)
% (Case A) y' = [y(i+1) - y(i-1)]/(2dx) error => (dx^2)
%
%
% setup matrix eqns. (treat boundary terms as special cases)
a = zeros(N,N); b = zeros(N,1);
for i = 2:N-1
a(i,i-1) = 1-3*x(i)*dx/2;
a(i,i) = -2+7*dx2;
a(i,i+1) = 1+3*x(i)*dx/2;
b(i) = dx2*cos(2*x(i));
end
% left boundary
a(1,1) = -2+7*dx2;
a(1,2) = 1+3*x(1)*dx/2;
b(1) = dx2*cos(2*x(1))-(1-3*x(1)*dx/2)*yxo;
% right boundary
a(N,N-1) = 1-3*x(N)*dx/2;
a(N,N) = -2+7*dx2;
b(N) = dx2*cos(2*x(N))-(1+3*x(N)*dx/2)*yxf;
% solve system of eqns
y = a\b;
% add boundary points to solution for plotting
za = [yxo y' yxf]; xa = [xo x xf];
%
%
% Approximations: y'' = [y(i-1) - 2y(i) + y(i+1)]/(dx^2) error => (dx^2)
% (Case B) y' = [y(i+1) - y(i)]/(dx) error => (dx)
%
%
% setup matrix eqns. (treat boundary terms as special cases)
a = zeros(N,N); b = zeros(N,1);
for i = 2:N-1
a(i,i-1) = 1;
a(i,i) = -(2+3*x(i)*dx-7*dx2);
a(i,i+1) = 1+3*x(i)*dx;
b(i) = dx2*cos(2*x(i));
end
% left boundary
a(1,1) = -(2+3*x(1)*dx-7*dx2);
a(1,2) = 1+3*x(1)*dx;
b(1) = dx2*cos(2*x(1))-yxo;
% right boundary
a(N,N-1) = 1;
a(N,N) = -(2+3*x(N)*dx-7*dx2);
b(N) = dx2*cos(2*x(N))-(1+3*x(N)*dx)*yxf;
% solve system of eqns
y = a\b;
% add boundary points to solution for plotting
zb = [yxo y' yxf]; xb = [xo x xf];
%
% plot results for both cases
t = '220+n';
subplot(eval(t)),plot(xa,za,'r-',xb,zb,'g--','LineWidth',2)
axis([0 3.2 -2 1]);
title(['FD: FD Method (',num2str(N),' pts)'])
xlabel('x values'),ylabel('y values'),grid
legend('Case A','Case B')
%
end
%
% end of demo