Plotting a temperature graphs of a heat equation of a rod

The values of c, L and deltat are choosen by myself. The equation is attached in the picture and this my code. I can't seem to get the curve correctly for the temperature of the rod. Is there anything wrong with my code?

syms n;

c = 4;

L = 6;

deltat = 0.2;

x = [0:6];

for k = 1:5

t = (deltat*k);

u = (200/pi)*symsum((((-1)^(n+1))/n)*(sin((n*x*pi)/L))*(exp(-(((c*n*pi)/L)^2)*t)),n,1,1000);

plot(x,u)

hold on

end

Answer 1

[Star Strider]

The Symbolic Math Toolbox isn’t the best option for this. It’s possible to vectorise the sin*exp term, and is only necessary to sum over ‘n’ in a loop:

c = 4;
L = 6;
deltat = 0.2;
k = 0:4;
t = k*deltat;
x = linspace(0,6,50);
xp = c*pi/L^2;
[T,X] = meshgrid(t,x);
u = zeros(size(T));
for n = 1:1000
    u = u + ((-1)^(n+1)/n) * (sin(n*pi*X/L).*exp(-n*xp*T));
end

figure(1)
surf(T, X, u)
grid on
xlabel('\itt\rm')
ylabel('\itx\rm')
zlabel('\itu(x,t)\rm')
view(125, 35)

Experiment to get the result you want.

[Benjamin Chee]

Thank you for your response. I really appreciate it. However, I have want question. There is a 200/pi inside the question. Can i include 200/pi into the equation ((-1)^(n+1)/n) * (sin(n*pi*X/L).*exp(-n*xp*T))?

[Star Strider]

My apologies. By the time I wrote the code for the summation, I forgot the 200/pi.

After the loop, multiply ‘u’ by 200/pi.

The loop (otherwise unchanged) will then be:

u = zeros(size(T));
for n = 1:1000
    u = u + ((-1)^(n+1)/n) * (sin(n*pi*X/L).*exp(-n*xp*T));
end
u = u*200/pi;

That multiplies the ‘u’ by 200/pi, as required.

[Benjamin Chee]

Thank you for the response. I actually tried putting 200/pi inside the loop equation. It also shows the similar result as the one placing 200/pi after the loop. Thank you for your help. I really appreciate it :)

[Star Strider]

As always, my pleasure.

If my Answer helped you solve your problem, please Accept it!

[Benjamin Chee]

Done ^^