Computing the Triple Summation of the Heat Equation

Hello, I have a general solution to my heat equation as:

I have tried implementing this in MATLAB but I get inconsistent results. Do you think the issue is about my code or just the triple summation does not yield realistic outputs. Here's my code.

% Define the bounds for the summation
n = 50; % Upper limit for the first index
m = 50; % Upper limit for the second index
p = 50; % Upper limit for the third index
to = 350;
tc = 80;
Lx = 10;
Ly = 10;
Lz = 10;
k = 1;

x=5;
y=5;
z=5;

% Initialize the sum variable
sum_val = 0;

% Perform the triple summation
for a = 1:n
    alpha = (2*a-1)*pi/Lx;
    for b = 1:m
        beta = (2*b-1)*pi/Ly;
        for c = 1:p
            gamma = (2*c-1)*pi/Lz;
            lambda = (alpha^2 + beta^2 + gamma^2);
            % Compute the term to be added to the sum
            term = (sin(alpha*x) * sin(beta*y) * sin(gamma*z) * exp(-lambda*k*0)) / ...
                ((2*a-1) + (2*b-1) + (2*c-1));

            
            % Add the term to the sum
            sum_val = sum_val + term;
        end
    end
end

u = to + sum_val*(64*(to-tc)/pi^3);

% Display the result
disp(['The triple summation value is: ' num2str(u)]);

[VBBV]

Can you tell how lambda is defined ? Also, the equation is computed at time t = 0 only. Consider heat equation that changes with time, to include the time loop along with spatial dimension.

[uk]

Each term on the right hand side also correspond to alpha, beta and gamma respectively. I set t=0 to initially see whether the code structure has flaws or not. t being set to 0 should not be a problem as of now I think, since the output does not correspond to initial conditions. Moreover, the output seems to decrease when I set a higher t value. So I guess the issue is something else?

[Torsten]

I wonder where rho and cp are hidden in your function U. Or does k equal alpha ?

[uk]

Yes, apologies for the confusion. K is indeed alpha.

[Torsten]

Seems you analytical solution as used is incorrect. Where did you find it ?

% Define the bounds for the summation
n = 50; % Upper limit for the first index
m = 50; % Upper limit for the second index
p = 50; % Upper limit for the third index
to = 350;
tc = 80;
Lx = 10;
Ly = 10;
Lz = 10;
k = 1;

x=5;
y=5;
z=5;

T = 0:1:100;
u = zeros(size(T));
% Perform the triple summation
for it = 1:numel(T)
  t = T(it); 
  % Initialize the sum variable
  sum_val = 0;
  for a = 1:n
    alpha = (2*a-1)*pi/Lx;
    for b = 1:m
      beta = (2*b-1)*pi/Ly;
      for c = 1:p
        gamma = (2*c-1)*pi/Lz;
        lambda = (alpha^2 + beta^2 + gamma^2);
        % Compute the term to be added to the sum
        term = (sin(alpha*x) * sin(beta*y) * sin(gamma*z) * exp(-lambda*k*t)) / ...
                ((2*a-1) * (2*b-1) * (2*c-1));

            
        % Add the term to the sum
        sum_val = sum_val + term;
      end
    end
  end
  u(it) = to + sum_val*(64*(tc-to)/pi^3);
end
plot(T,u)

[uk]

In a book called Introduction to Partial Differential Equations by Russel Herman. Here's the link if you want to check it out. The section starts at page 185. You can also check against the diagrams produced by the author.

[uk]

It seems that the issue was due to the fact that we added terms in the denominator but they need to be multiplied. The results are still unrealistic though.

[Torsten]

Seems to work now (see above).

Additionally to the multiplication in the denominator the line

u(it) = to + sum_val*(64*(to-tc)/pi^3);

had to be changed to

u(it) = to + sum_val*(64*(tc-to)/pi^3);

But you should check against a numerical code to see if there might be other errors.