How can I change my code so that the c matrix is not overwritten
2 次查看(过去 30 天)
显示 更早的评论
This is my code modelling advection diffusion to obtain circular diffusion graphs, i get the error Index in position 3 exceeds array bounds (must not exceed 1). I understand I get this because for line [C,h]=contourf(X,Y,Z,300); I am changing the matrix into 2d, and therfore the 3rd dimension is non existent. How can I change the code so that I dont have to overwrite the matrix or that the graphs plot at the same time before the new matrix is assigned.
clf
clear
clc
% Inputs
D=1.04e-6; % Thermal diffusivity, m2/s, eg 0.01
r_in=0;
r_out=0.05;
m=35; %number of divided sections between r_out and r_in
delta_r=(r_out-r_in)/m; %section length, m
nr=m+1; %total no. of radial points
nphi=36;%no. of angle steps
delta_phi=2*pi/nphi; %angle step
t=40000; % total time, s
u=delta_r/t; % veloctiy in x direction, m/s
nt=450000; % total no. of time steps
delta_t=t/nt; % timestep, s
Crin=0; % Concentration at r_in, ppm
Crout=1000; % Concentration at r_out, ppm
Cin=0; %initial concentration at r_in<r<r_out, ppm
Cmax=max([Crin,Crout]);
% Solution
for i=1:nr
r(i)=r_in+(i-1)*delta_r;
end
r;
disp(u)
phideg=[0:360/nphi:360+360/nphi];
phi=phideg.*pi/180;
c=u*delta_t/delta_r; % Courant Number
d=D*delta_t/delta_r^2; % diffusion number
d1=D*delta_t/(2*delta_r);
Pe=u*r/D; % Peclet Number
%fprintf('\n')
if d<=0.5
fprintf('solution stable\nd=%10.7f\n', d)
else
fprintf('solution unstable\nd=%10.7f\n', d)
end
fprintf('\n')
fprintf('\n')
if c^2<=2*d
fprintf('solution stable\n c^2 (%4.2f)<=2*d(%4.2f)\n', c^2, 2*d)
else
fprintf('solution unstable\n c^2 (%4.2f)>(%4.2f)\n', c^2, 2*d)
end
% Solution
for i=1:nr
r(i)=r_in+(i-1)*delta_r;
end
% Creating initial boundary conditions
C=zeros(nr,nphi+2,nt);
for k = 1:nt+1
for j = 1:nphi+2
for i=1:nr
if(i==1)
C(i,j,k)=Crin;
elseif(i==nr)
C(i,j,k)=Crout;
else
C(i,j,k)=Cin;
end
end
end
end
C;
% Creating C matrix
for k=1:nt
for j=1:nphi+2
for i= 2:nr-1
C(i,j,k+1)= C(i,j,k) +(c/2)*(C(i+1,j,k)-C(i-1,j,k)) + d*(C(i+1,j,k) - 2*C(i,j,k) + C(i-1,j,k))+(d1/r(i))*(C(i+1,j,k)-C(i-1,j,k));
end
end
end
C(:,:,1);
C(:,:,nt+1);
C;
C(:,1,1);
C(:,1,nt+1);
C(:,1,:);
%Initial Temperature profile
C1=C(:,:,1);
subplot(2,2,1)
r=r;
[Phi,R]=meshgrid(phi,r);
Z=C1; %C1.*Phi./Phi;
[X,Y,Z]=pol2cart(Phi,R,Z);
%figure(1)
[C,h]=contourf(X,Y,Z,300);
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel(['r in m'])
%final C
C2=C(:,:,nt+1); % <-------------- Where the error appears
subplot(2,2,3)
r=r;
[Phi,R]=meshgrid(phi,r);
Z=C2.*Phi./Phi;
[X,Y,Z]=pol2cart(Phi,R,Z);
%figure(1)
[C,h]=contourf(X,Y,Z,300);
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel(['r in m'])
1 个评论
回答(1 个)
Maneet Kaur Bagga
2024-5-8
Hi,
As per my understanding the error is encountered when you are trying to visualize the results of a 3D matrix "C" as "2D" plots at different time steps. The error is because you're trying to access a slice of the matrix "C" that does not exist in a 2D context after transforming it via "pol2cart". The transformation and contour plotting overwrite the "C" variable.
To resolve this, you should avoid overwriting the variable C with the contour plot handle. In MATLAB, [C,h]=contourf(...) assigns the contour matrix to C, which is also the name of your data matrix. This overwrites your original C matrix when you meant to only create a contour plot.
The possible workaround for this is:
Change the variable names for contour plots: Use a different variable name for the contour plot handles to avoid overwriting the "C" matrix. For example, use [Ch,h] instead of [C,h] when calling contourf.
Please refer to the code snippet below for reference:
% Initial Temperature profile
C1 = C(:,:,1);
subplot(2,2,1)
[Phi,R] = meshgrid(phi,r);
Z = C1; % No need to modify Z here
[X,Y,Z] = pol2cart(Phi,R,Z);
[Ch,h] = contourf(X,Y,Z,300); % Changed C to Ch to avoid overwriting
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel('r in m')
% Final C
C2 = C(:,:,nt+1);
subplot(2,2,3)
[Phi,R] = meshgrid(phi,r);
Z = C2; % Directly use C2, assuming you want to plot it as is
[X,Y,Z] = pol2cart(Phi,R,Z);
[Ch,h] = contourf(X,Y,Z,300); % Changed C to Ch to avoid overwriting
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Final Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel('r in m')
Hope this helps!
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Contour Plots 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!