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:   -(x2 y'(x))' + x3/(2+x) y(x) = x2/(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
>> xh = a+h : h : b-h;                % interior nodes
>> xp = (a + 0.5*h) : h : (b-0.5*h);  % "half-nodes"
>> hq = h*h * xh .^ 3 ./ (2 + xh);    % node values of q, multiplied by h*h
>> p = xp .^ 2;                       % "half-node" values of p
>> n = length(xh);                    % the number of unknowns
>> 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*h * xh .^ 2 ./ (3 - xh);  % node values of f, multiplied by h*h
>> 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,'ko')          % a graph of the solution

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

Self-adjoint formulation:
-(e2x y'(x))' - e2x y(x) = -ex

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
>> xh = a+h : h : b-h;                % interior nodes
>> xp = (a + 0.5*h) : h : (b-0.5*h);  % "half-nodes"
>> hq = -h*h *exp(2*xh);              % node values of q, multiplied by h*h
>> p = exp(2*xp);                     % "half-node" values of p
>> n = length(xh);                    % the number of unknowns
>> 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*h *exp(xh) ;        % node values of f, multiplied by h*h
>> 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,'ko')         % a graph of the solution

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

graphically - red asterisks represent the exact solution at the nodes

>> yy = (2 - 2*x + 0.5*x.*x).*exp(-x);
>> hold on
>> plot(x,yy,'r*')

numerically - compute the error

>> max(abs(y-yy'))
  ans =  5.7933e-04