error estimate using finite difference method

% Program to implement the 
% Backward in time-Central in space
% Method to solve heat conduction equation
        
clear all
close all

format long
clc
% Input variables
%alpha = 1;
N = 15; % Number of subintervals

error = zeros(5,1);
order = zeros(4,1);
h = zeros(5,1);
% Calculate step size
for p=1:5
    M = 1000; % Time step size
    a = 0; b = 1; % Domain
    t0 = 0; tf = 1;
    h(p) = (b-a)/N;
    k = (tf-t0)/M;
    % Unconditionally stable
    mu = k/(h(p)^2);
for i=1:N+1
        x(i) = a + (i-1)*h(p);
f= sin(pi*x) + sin(2*pi*x);
end
u = zeros(N+1,M);
% Find the approximate solution at each time step
for n = 1:M
t = n*k; % current time
% boundary condition at left side
gl = sin(pi*a)*exp(-pi*pi*t)+sin(2*pi*a)*exp(-4*pi*pi*t);
% boundary condition at right side
gr = sin(pi*b)*exp(-pi*pi*t)+sin(2*pi*b)*exp(-4*pi*pi*t);
if n==1 % first time step
for j=2:N % interior nodes
u(j,n) = f(j) + mu*(f(j+1)-2*f(j)+f(j-1));
end
u(1,n) = gl; % the left-end point
u(N+1,n) = gr; % the right-end point
else
for j=2:N % interior nodes
u(j,n)=u(j,n-1)+mu*(u(j+1,n-1)-2*u(j,n-1)+u(j-1,n-1));
end
u(1,n) = gl; % the left-end point
u(N+1,n) = gr; % the right-end point
end
end
%u(q) = u(j,n);
N = N*2;
%for j = 1:N-1
 error(p) = norm((u(j+1,n)-u(j,n)),2);
   % N = N*2;
%end
    for j=1:p-1
        order(j) = log(error(j)/error(j+1))/log(2);
    end
end

I need help in error statement. I want to estimate error from previous value of h to current value of h.

[Torsten]

The line

error(p) = norm((u(j+1,n)-u(j,n)),2);

doesn't make sense since you take the norm of a single value, namely u(j+1,n)-u(j,n).

Which error norm of the discrete approximations do you try to compute ?

[Kumar]

I have defined h(p) = (b-a)/N; and later i have taken N = 2*N which gives h = [0.0667, 0.0333, 0.0167, 0.0083, 0.0042]. I want to approximate error at next values of h to current value of h(i.e, for e.g, error = (u for h = 0.0333) - (u for h = 0.0667) and so on).

[Torsten]

I understood this, but you didn't specify how you define

(u for h = 0.0333) - (u for h = 0.0667)

Only every second value from (u for h = 0.0333) exists for (u for h = 0.0667).

Then say you subtract the vectors for u at points that both calculations share, you want to take the 2-norm of this difference vector to define the error ?

[Kumar]

yes, i want to take 2 norm to estimate the error.

[Torsten]

I plotted the solution curves for h=1/15 for start and end time of your simulation.

For which output time(s) do you want to compute the differences in the solution for different values of h ?

% Program to implement the 
% Backward in time-Central in space
% Method to solve heat conduction equation
        
clear all
close all

format long
clc
% Input variables
%alpha = 1;
N = 15; % Number of subintervals

error = zeros(5,1);
order = zeros(4,1);
h = zeros(5,1);
% Calculate step size
for p=1:1
    M = 1000; % Time step size
    a = 0; b = 1; % Domain
    t0 = 0; tf = 1;
    h(p) = (b-a)/N;
    k = (tf-t0)/M;
    % Unconditionally stable
    mu = k/(h(p)^2);
    for i=1:N+1
        x(i) = a + (i-1)*h(p);
        f= sin(pi*x) + sin(2*pi*x);
    end
    u = zeros(N+1,M);
    % Find the approximate solution at each time step
    for n = 1:M
       t = n*k; % current time
       % boundary condition at left side
       gl = sin(pi*a)*exp(-pi*pi*t)+sin(2*pi*a)*exp(-4*pi*pi*t);
       % boundary condition at right side
       gr = sin(pi*b)*exp(-pi*pi*t)+sin(2*pi*b)*exp(-4*pi*pi*t);
       if n==1 % first time step
          for j=2:N % interior nodes
             u(j,n) = f(j) + mu*(f(j+1)-2*f(j)+f(j-1));
          end
          u(1,n) = gl; % the left-end point
          u(N+1,n) = gr; % the right-end point
       else
          for j=2:N % interior nodes
             u(j,n)=u(j,n-1)+mu*(u(j+1,n-1)-2*u(j,n-1)+u(j-1,n-1));
          end
          u(1,n) = gl; % the left-end point
          u(N+1,n) = gr; % the right-end point
      end
   end
   %u(q) = u(j,n);
    N = N*2;
   %for j = 1:N-1
    error(p) = norm((u(j+1,n)-u(j,n)),2);
   % N = N*2;
   %end
   for j=1:p-1
      order(j) = log(error(j)/error(j+1))/log(2);
   end
end
plot(x,[u(:,1),u(:,end)])