How can I change my code so that the c matrix is not overwritten
2 views (last 30 days)
Show older comments
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 Comment
Answers (1)
Maneet Kaur Bagga
on 8 May 2024
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 Comments
See Also
Categories
Find more on Contour Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!