1

I have a general solution to my heat equation as: triple summation 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)]);

The code outputs a value of 403.4368 for x,y,z = 5. It should not be this high. I was expecting a result of 80 at most.

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
Orkiel
  • 31
  • 2
  • The equation says to sum up to infinity, you sum up only to 50. I haven’t looked at the function, I hope it is very close to zero after 50. – Cris Luengo Jul 02 '23 at 16:44

1 Answers1

2

I looked over the code, and noticed a typo - you have a sum in the denominator, instead of a product. So, instead of

term = (sin(alpha*x) * sin(beta*y) * sin(gamma*z) * exp(-lambda*k*0)) / ...
            ((2*a-1) + (2*b-1) + (2*c-1));

according to your equation it should be

term = (sin(alpha*x) * sin(beta*y) * sin(gamma*z) * exp(-lambda*k*0)) / ...
            ((2*a-1) * (2*b-1) * (2*c-1));
Brionius
  • 13,858
  • 3
  • 38
  • 49