Partial differential equation, error estimation

Hi there,

Im busy with a nummerical method for approximating a partial differential equation. The PDE is the temperature in a bar. The partial differential equation that describes the scaled temperature in the bar is \begin{center}

delta T / delta t = a^2 delta^2T / dellta x^2

In which a^2 is the thermal diffusitivity, which has for our bar a value of a^2 = 0.25. The bar covers the interval x in range [0, 1.5]. All the parameters are well fixed in the matlab script. We obtain the following boundary conditions:

T(t,0) = 1 + beta sin(t),

T(t,1.5) = 1

The intitial condition is given by

T(0,x) = 1

First I used semi discretization, then trapezoidal method to estimate the temperature in the bar after 10 seconds. From there i wanna calculate the global truncation error. however that went wrong. Delta t is given by:

dt = (5*dx2)/(a.^2);

N = ceil(T/dt2)

dt = T/N;

Is there someone who can help me?

clear all; close all; clc;

% setting parameters
a = sqrt(0.25);
ti = 0;
tf = 10;
RB = 1;
beta = 0.05;
xi=0;
xf=1.5;

% setting initial parameters (those will change in the for loop)

k = [0; 1; 2; 3; 4; 5];
for j = 1:length(k)

dx1 = (1.5)/(10*2^k(j));
dt1 = (5*dx1)/(a.^2);
T = 10;
N = ceil(T/dt1);
dt1 = T/N;


x1 = [xi+dx1:dx1:xf-dx1];
n = round((xf-xi)/(dx1))-1;
K1=-2*eye(length(x1))*((a.^2)/(dx1^2)); 
for i= 1: length(x1)-1
    K1(i,i+1)= ((a.^2)/(dx1^2)) ;
    K1(i+1,i)= ((a.^2)/(dx1^2)) ;
end

t1 = linspace(dt1,tf, (tf/dt1));
w1 = zeros(length(x1),length(t1));
r1 = zeros(length(x1),1);
rjplus1 = zeros(length(x1),1);

for n1 = 1:length(t1)
    LB1 = 1+beta*sin(t1(n1));
    r1(1,1) = 1+beta*sin(t1(n1)-dt1); r1(end) = 1;
    r1 = r1*((a.^2)/(dx1.^2));
    
    rjplus1(1,1) = 1+beta*sin(t1(n1)); rjplus1(end) = 1;
    rjplus1 = rjplus1*((a.^2)/(dx1.^2));
    
    A1 = eye(length(x1))-(dt1/2)*K1;
    B1 = w1(:,n1)+(dt1/2)*(K1*w1(:,n1)+r1(:)+rjplus1(:));
    w1(:,n1+1) = A1\B1;
end

T10_1 = [LB1; w1(:,end); RB]; % Redefine temperature with boundaries added
Xextended1 = [0; x1(:); RB];
AvgT1 = trapz(Xextended1,T10_1); % avg scaled temperature from trapezoidal rule

dx2 = 2*dx1; % defining the second data set with dx = 2*dx
dt2 = (5*dx2)/(a.^2);
T = 10;
N = ceil(T/dt2);
dt2 = T/N;

x2 = [xi+dx2:dx2:xf-dx2];

K2=-2*eye(length(x2))*((a.^2)/(dx2^2)); 
for i= 1: length(x2)-1
    K2(i,i+1)= ((a.^2)/(dx2^2)) ;
    K2(i+1,i)= ((a.^2)/(dx2^2)) ;
end

t2 = linspace(dt2,tf, (tf/dt2));
w2 = zeros(length(x2),length(t2));
r2 = zeros(length(x2),1);
rjplus2 = zeros(length(x2),1);

for n2 = 1:length(t2)
    LB2 = 1+beta*sin(t2(n2));
    r2(1,1) = 1+beta*sin(t2(n2)-dt2); r2(end) = 1;
    r2 = r2*((a.^2)/(dx2.^2));
    
    rjplus2(1,1) = 1+beta*sin(t2(n2)); rjplus2(end) = 1;
    rjplus2 = rjplus2*((a.^2)/(dx2.^2));
    
    A2 = eye(length(x2))-(dt2/2)*K2;
    B2 = w2(:,n2)+(dt2/2)*(K2*w2(:,n2)+r2(:)+rjplus2(:));
    w2(:,n2+1) = A2\B2;
end

T10_2 = [LB2; w2(:,end); RB]; % Redefine temperature with boundaries added
Xextended2 = [0; x2(:); RB];
AvgT2 = trapz(Xextended2, T10_2); % avg scaled temperature from trapezoidal rule

xv1 = round(length(x1)/2); % selecting the value when x = 0.75 in order to calculate the error
xv2 = round(length(x2)/2); 

p = 2; %because it is second order
error(j) = abs(T10_1(xv1)-T10_2(xv2))/((2^p)-1);


end

N = zeros(length(k),1);
for i = 1:length(k)

dx(i) = (1.5)/(10*2^k(i));
dt1(i) = (5*dx(i))/(a.^2);
T(i) = 10;
N(i) = ceil(T(i)/dt1(i));
dt(i) = T(i)/N(i);
end

% creating a table with k, dt, dx and the error
dx = dx'; dt = dt'; error = error';
table = table(k, dt, dx, error)
dxc = dx(6);
kxc = k(6);
txc = dt(6);
formatSpec = 'The absolute value of the global truncation error in T(t) at location x = 0.75 with beta = 0.05 is smaller than 10^-5 when k = %1f, dx is %4.3f and dt is %4.3f\n' ;
fprintf(formatSpec,kxc,dxc, txc)

% plotting the errors in a nice graph
limit = ones(length(dx))*10^(-5); 
loglog(dx,error,'-bp')
hold on
loglog(dx,limit)
title('Global truncation error per different values dx and so dt and k')
xlabel('dx')
ylabel('global truncation error')
legend('error','10E-5')

[J. Alex Lee]

What exactly do you mean "global truncation error", and what exactly is "going wrong"?

What are T10_1 and T10_2?

Can you lay out your pencil-and-paper math reasoning for how to calculate what you are trying to calculate?