Finite difference method in 1D

Second order linear ODE (selfadjoint formulation)

Second order linear ODE in selfadjoint formulation:

-(p(x) y'(x))' + q(x) y(x) = f(x) in (a,b), y(a) = y0, y(b) = yn .

Problem:   -(x² y'(x))' + x³/(2+x) y(x) = x²/(3-x) ,    y(-5) = -2 ,   y(-3) = 0 .

Solving the problem in Matlab by finite difference method with step-size h = 0.4:

 % assembling a system matrix A
 h = 0.4;                           % step-size
 a = -5; b = -3;                    % interval of solution
 n = round((b-a)/h)-1;              % the number of unknowns
 xh = linspace(a+h, b-h, n);        % interior nodes
 xp = (a + 0.5*h) : h : b;          % "half-nodes"
 hq = h^2 * xh .^ 3 ./ (2 + xh);    % node values of q, multiplied by h^2
 p = xp .^ 2;                       % "half-node" values of p
 D = hq + p(1:n) + p(2:n+1);        % the main diagonal
 V = - p(2:n);                      % the second diagonal
 A = diag(V,-1) + diag(D) + diag(V,1);  % system matrix

 % the right hand side
 y0 = -2; yn = 0;                % Dirichlet boundary values
 F = h^2 * xh .^ 2 ./ (3 - xh);  % node values of f, multiplied by h^2
 F(1)=F(1) + p(1)*y0;            % incorporating the left boundary value
 F(n)=F(n) + p(n+1)*yn;          % incorporating the right boundary value

 % computation of the solution and graph plotting
 U = A\F';                   % solution of the system
 y = [y0; U ; yn];           % extension to boundary values
 x = [a,xh,b];               % extension to end-points of the interval
 plot(x,y,x,y,'ro')          % graph of the numerical solution

Problem - dumped oscillations:   y''(x) + 2 y'(x) + y(x) = e^-x,    y(0) = 2 ,   y(2) = 0 .

Self-adjoint formulation:
-(e^2x y'(x))' - e^2x y(x) = -e^x

Solving the problem in Matlab by finite difference method with step-size h = 0.2:

 %  assembling a system matrix A
 h = 0.2;                           % step-size
 a = 0; b = 2;                      % interval of solution
 n = round((b-a)/h)-1;              % the number of unknowns
 xh = linspace(a+h, b-h, n);        % interior nodes
 xp = (a + 0.5*h) : h : b;          % "half-nodes"
 hq = -h^2 *exp(2*xh);              % node values of q, multiplied by h^2
 p = exp(2*xp);                     % "half-node" values of p
 D = hq + p(1:n) + p(2:n+1);        % the main diagonal
 V = - p(2:n);                      % the second diagonal
 A = diag(V,-1) + diag(D) + diag(V,1);  % system matrix

 % the right hand side
 y0 = 2; yn = 0;            % Dirichlet boundary values
 F = -h^2 *exp(xh) ;        % node values of f, multiplied by h^2
 F(1) = F(1) + p(1)*y0;     % incorporating the left boundary value
 F(n) = F(n) + p(n+1)*yn;   % incorporating the right boundary value

 % computation of the solution and graph plotting
 U = A\F';                  % solution of the system
 y = [y0; U ; yn];          % extension to boundary values
 x = [a,xh,b];              % extension to end-points of the interval
 plot(x,y,'ro')             % graph of the numerical solution

compare the results with the exat solution of y = (2 - 2x + 0.5 x²) e^-x

graphically - the black curve

 xx = linspace(a,b,30);
 yy = (2 - 2*xx + 0.5*xx.^2).*exp(-xx);
 hold on
 plot(xx,yy,'k')

compare the results numerically - compute the error

 yy = (2 - 2*x + 0.5*x.^2).*exp(-x);
 max_err = max(abs(y-yy'))