clear
tic
L=1;
T=0.2;
nust=2000;
dt=T/nust;
n=40;
dx=L/n;
r=1;
omega=10:10:5000;%Store Range of Frequencies for Simulation
u=zeros(n+1,nust+1);%Initialize the grid (Space along rows/dim1& time along cols/dim2)
% Bar Boundary Condition
time=zeros(1,nust+1);%Time Elapsed Vector
for k=1:nust+1
u(1,k)=0;
u(n+1,k)=0;
time(k)=(k-1)*dt;
end
% Bar Initial Condition
x=zeros(n+1,1);%Length Vector
fori=1:n+1
x(i)=(i-1)*dx;
u(i,1)=sin(pi*x(i));
end
%constructing the Matrix in the right side
ar(1:n-2)=r;
br(1:n-1)=2-2*r;
cr(1:n-2)=r;
Mr=diag(br,0)+diag(ar,-1)+diag(cr,1);
%constructing the Matrix in the left side
al(1:n-2)=-r;
bl(1:n-1)=2+2*r;
cl(1:n-2)=-r;
Ml=diag(bl,0)+diag(al,-1)+diag(cl,1);
%aa=Ml\Mr;
% Crank Nicolson Implementation
u_sol=zeros(length(u(:,1)),length(u(1,:)),length(omega));
forll=1:length(omega)
for k=2: nust
C=0.5*dt*(cos(omega(ll)*k*dt)+cos(omega(ll)*(k+1)*dt))*ones(length(2:n),1);
uu=u(2:n,k-1);
u(2:n,k)=inv(Ml)*((Mr*uu)+C);
end
u_sol(:,:,ll)=u;
clear uu;
end
figure;
subplot(2,1,1); mesh(x,time,u_sol(:,:,1)')
title(['Solution for omega=',num2str(omega(1))]);
subplot(2,1,2); mesh(x,time,u_sol(:,:,end)')
title(['Solution for omega=',num2str(omega(end))]);
xlabel('x-axis');
ylabel('Temperature');
toc
The code you have provided me at that time was giving the out of the whole process( as a plot ) which isgood but we cannot specify a point in the rod to find the temperature , for example if we need to know the temperature at the distance 0.5 ( we know the total length is 1) there is no way to get it . so what I need is the following :
We have a rod with a length of 1 and we want to know the temperature at the center of the rod ( we may change the point we want to know ) and the output must be a number not a plot , I mean I want the MATLAB to tell me the temperature at the specified point.