Matlab code for the Formula

Following is the Analytical solution of a Heat Diffusion formula where x is length of Rod from 0 to 1 meter with an increment of 0.1. Also t is the time from 0 to 10 minute with an increment of 1 min. Value of r is 0.1. The result should almost match with the following table
. The result in the table was obtained using the Numerical Method (Explicit Method) of the problem.
I've tried the following code:
close all clear clc
x = 0:0.1:1; t = 0:1:10; r = 0.1; [X,T] = meshgrid(x,t) S = 0;
for n=1:1e6 S = S + sin(pi*n/2)*sin(n*pi*X).*exp(-n*n*r*r*t)/(n*n); end U = (8/pi*pi)*S
What is wrong with this code?

2 个评论

What have you tried? If nothing, why not? It is by making an effort that you will learn, not by being given the solution to your homework on a platter. So show what you tried. If you do, you will have abetter chance of getting some help here.
close all
clear
clc
x = 0:0.1:1;
t = 0:1:10;
r = 0.1;
[X,T] = meshgrid(x,t)
S = 0;
for n=1:1e6
S = S + sin(pi*n/2)*sin(n*pi*X).*exp(-n*n*r*r*t)/(n*n);
end
U = (8/pi*pi)*S

请先登录,再进行评论。

回答(2 个)

One way is a more direct, though not vectorized, approach of using for loops over t and x:
x = 0:0.1:1;
t = 0:1:10;
r = 0.1;
U = zeros(length(t), length(x));
for kt = 1 : length(t)
for kx = 1 : length(x)
s = 0;
for n=1:1e6
s = s + sin(pi*n/2).*sin(n*pi*x(kx)).*exp(-n*n*r*r*t(kt))/(n*n);
end
U(kt, kx) = 8 * s / pi^2;
end
end
U
Hey Rasel, your code was almost right. You made just a few errors, instead of initializing the sum as you did, you should have just initialized U as a matrix. This is your code with a few changes.
clear variables
close all
x = 0:0.1:1;
t = 0:10;
r = 0.1;
N = length(x);
M = length(t);
U = zeros(N,M); % This is what you should have done to start with.
[x,t] = meshgrid(x,t);
for n = 1: 100
U = U + ((1/n^2)*sin(0.5*pi*n)*sin(pi*n*x).*exp(-n^2*r^2*t))';
end
U = (8/pi.^2)*U;
disp(U')
figure
surf(x,t,U')
grid

类别

帮助中心File Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by