I do not know what's wrong with my code

1 view (last 30 days)
Nathan Bueno
Nathan Bueno on 12 Mar 2019
Answered: Steven Lord on 12 Mar 2019
%1-D Heat equation
%example 1 at page 782
%dT/dt=c.d^2T/dx^2
%T(x,t)=temperature along the rod
%by finite difference method
%(T(x,t+dt)-T(x,t))/dt = (T(x+dx,t)-2T(x,t)+T(x-dx,t))/dx^2
%solve for T(x,t+dt) by iteration
%heat constant
close all;
%clear all;
c=1;
L=1; %length of the rod
N=5; %# of elements
%dicretize Xspace
x_vec=0:0.2:1;
dx=x_vec(2)-x_vec(1);
%discretize time
dt=0.5/50;
t_vec=0:dt:0.5;
%temperature matrix
T_mat=zeros(length(x_vec),length(t_vec));
T_mat1=zeros(length(x_vec),length(t_vec));
%boundary conditions
T_mat(1,:)=0;
T_mat(end,:)=0;
%initial conditions
T_mat(:,1)= 1-cos(pi*x_vec(idx))(400/npi)*sin*(npi.*x_vec(idx));
%T_mat(:,1)= sin*((pi.*x_vec);
%[tt,xx]=meshgrid(t_vec,x_vec);
subplot(2,1,1);
mesh(xx,tt,T_mat);
lamda=(c*dt/(dx^2));
for tdx=1:length(t_vec)-1
for idx=1:length(x_vec)
if idx==1
T_mat(idx,tdx+1)=T_mat(idx+1,tdx);
T_mat1(idx,tdx)= 1-cos(pi*x_vec(idx))(400/npi)sin(npix_vec(idx))*exp(-1*(pi^2)*t_vec(tdx)))); %modified code, 2nd line of 'if'
elseif idx==length(x_vec)
T_mat(idx,tdx+1)=T_mat(idx-1,tdx);
T_mat1(idx,tdx)= 1-cos(pi*x_vec(idx))(400/npi)sin*(npi*x_vec(idx))*exp(-1*(pi^2)*t_vec(tdx)); %modified code, 2nd line of 'elseif
else
T_mat(idx,tdx+1)=T_mat(idx,tdx)+lamda*((T_mat(idx+1,tdx)-2*T_mat(idx,tdx)+T_mat(idx-1,tdx)));
T_mat1(idx,tdx)=1-cos(pi*x_vec(idx))(400/npi)*sin(npix_vec(idx))*exp(-1*(pi^2)*t_vec(tdx)); %modified code, 2nd line of 'elseif)
end
end
end
%plot
subplot(2,1,1)
[tt,xx]=meshgrid(t_vec,x_vec);
mesh(xx,tt,T_mat);
title('Finite')
%plot
subplot(2,1,2)
[tt,xx]=meshgrid(t_vec, x_vec);
mesh(xx,tt,T_mat1);
title('Analytical');
figure
plot(x_vec,T_mat(:,1),x_vec,T_mat(:,11),x_vec,T_mat(:,21),x_vec,T_mat(:,31),x_vec,T_mat(:,41),x_vec,T_mat(:,51));
xlabel('Rod length (m)'); ylabel('Temperature (C)');
legend('Initially','At t=0.1sec','At t=0.2 secs','At t=0.3 secs',' At t=0.4 secs',' At t=0.5 secs');
figure
plot(x_vec,T_mat1(:,1),x_vec,T_mat1(:,11),x_vec,T_mat1(:,21),x_vec,T_mat1(:,31),x_vec,T_mat1(:,41),x_vec,T_mat1(:,51));
xlabel('Rod length (m)'); ylabel('Temperature (C)');
legend('Initially','At t=0.1sec','At t=0.2 secs','At t=0.3 secs',' At t=0.4 secs',' At t=0.5 secs');
title('Analytical')
xlabel('Rod Length |m'); ylabel('Temperature |C|');
legend('Initially','At t=0.1 sec','At t=0.0 sec','At t=0.3 sec','At t=0.4 sec','At t=0.5 sec');
  1 Comment
Raghunandan V
Raghunandan V on 12 Mar 2019
https://in.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

Sign in to comment.

Answers (1)

Steven Lord
Steven Lord on 12 Mar 2019
If you receive an error when running your code, post the full text (everything displayed in red) of the error message you receive, and indicate (if it's not obvious) the line in your code on which the error occurs.
If you don't receive an error but receive a different result than you expected, it's likely going to be difficult for us to help you since there aren't a lot of comments in the code to help us understand it. This looks like it may be the solution to a homework assignment (the "%example 1 at page 782" suggests that to me) and if that's true you may want to ask your professor and/or teaching assistant to help you. If it's not a homework assignment, or if your professor and/or TA are unavailable, step through your code using the debugging tools and identify where the code gives a different result than you expected then figure out why.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!