initial boundary condition for neumann conditions

clear;
clc;

format long;

% Thermal Diffusivity
alpha = 0.1;

% Length and Temperatures
L = 1;

% Number of Division for Space and Time
N_space = 21;
N_time = 11;

% Equispaced Grids for Space and Time
delta_x = 0.05;
delta_t = 0.01;

x = linspace(0, L, N_space)';
t = linspace(0, 0.5, N_time)';

d = alpha*delta_t/(delta_x)^2;

% initial conditions
U = x;

% boundary conditions
U(2) = U(1);
U(N_space - 1) = U(N_space);

% FTSC Explicit Scheme

c_0 = L;
tolerance = 10^(-12);
converged = 0;
m = 0;
result = 0;
while ~converged
    n = 2*m + 1;
    c_n = (2*L)/(n*pi)^2*((-1)^n - 1);
    T = c_n*exp(-0.5*alpha*((n*pi)/L)^2)*cos((n*pi*x)/L);
    result = result + T ;
    m = m + 1;
    if abs(T/result) < tolerance
        converged = 1;
    end
end
T = (c_0/2) + result;

subplot(2,1,1);
for k = 1:N_time - 1
    for i = 2:N_space - 1 % for loop indexes points inside domain only!
        V(i) = U(i) + d*(U(i + 1) - 2*U(i) + U(i - 1));
    end
    % impose boundary conditions
    V(1) = V(2);
    V(N_space) = V(N_space - 1);
    U = V'; % increment the solution in time
    if k == 3 - 1 || k == 5 - 1 || k == 7 - 1 || k == 9 - 1 || k == 11 - 1
        plot(x, U, 'DisplayName', ['t = ' num2str(t(k+1))]);
        hold on
        legend('show')
        title('FTCS Explicit: dx = 0.05, dt = 0.05');
        xlabel('L (ft)');
        ylabel('T (F)');
    end
    if k == 11 - 1
        abserr = abs(U - T);
        subplot(2,1,2);
        plot(x, abserr);
        title('Error Analysis of FTCS Explicit at t = 0.5 hr: dx = dt = 0.05')
        xlabel('L (ft)');
        ylabel('Error');
    end
end

I have a part in the code that I am curious about: my boundary conditions before my for loop. I am to code the numerical methods to solve the heat equation with specific initial conditions, and I was wondering what is the right way to code it.

I know the initial conditions for U(1) and U(2) before the for loop, but if we apply the first-order forward difference, then we must say that U(1) = U(2). However, knowing that coding is case-sensitive, the way to write your initial conditions will change the outcome. Should it be more accurate to say U(1) = U(2), or U(2) = U(1) before my for loop?

For example, U(1) = 0 and U(2) = 1 for my initial conditions. However, before my for loop, should I set it as

U(1) = 1 and U(2) = 1, or

U(1) = 0 and U(2) = 0?

Apparently for the FTCS explicit method code that I am showing here, if I set delta t = 0.01 instead, then my error analysis says that setting U(1) = 0, U(2) = 0 before the for loop gives a more accurate curve than setting U(1) = 1, U(2) = 1 before the for loop. Both situations end up with similar curves, but there is less error when U(1) = U(2) = 0 than U(2) = U(1) = 1 before the for loop.

I even have codes for FTCS implicit method and Crank-Nicolson method, and both of them is saying that there is less error when U(1) = U(2) = 0 than when U(2) = U(1) = 1 before the for loop.

Another reason why I am asking this question is because, in our for loop, we need to recalculate U(1) and U(2); however, our for loop calculates U(2) only, so the only way to give U(1) a value is, by using the first-order forward difference method, setting U(1) = U(2).

For example, supposed that U(2) = 0.3 after one computation. U(1) is now unknown. So the only way for U(1) to have a value is set U(1) = U(2) = 0.3.

As you can see, in our for loop, we have no choice but to code it as U(1) = U(2), but before the for loop, we have a choice to code it as U(1) = U(2) or U(2) = U(1). The only thing is which should we choose for prior to the for loop?

U(1) = U(2) (more consistent with the for loop), or

U(2) = U(1) (ends up with a more accurate curve)

[J. Alex Lee]

Potentially a weakness of your code is that your ought to be 2nd order accurate, but the way you apply your BC (zero flux) is only first order accurate.

To some extent your problem seems contrived because you area simultaneously wanting to assert an initial condition that looks like a striaght line (presumably), but that still satisfies zero-gradient at the ends - a contradiction that you can only tweak numerically.

I think the way you are evaluating "accuracy" is thus dubious.

Maybe it will be helpful to stop thinking of the BC's as "equating nodal values", and think more as that there is an additional relationship between the edge nodes that is not the same as your internal FDEs. For example, you DO have a choice on how to discretize your BC. Why not a 2nd order accurate 3-node derivative expression to relate the first 3 nodes rather than the first 2? How will your thinking change then?