How can I change my code so that the c matrix is not overwritten

1 次查看（过去 30 天）
Marjaan Chowdhury
Marjaan Chowdhury 2021-5-23
回答： Maneet Kaur Bagga 2024-5-8
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 个评论
Torsten
Torsten 2021-5-23
编辑：Torsten 2021-5-23
Why not naming the matrix in the call to contourf differently (e.g. M) ?
Why do you need the contour matrix at all ?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Maneet Kaur Bagga
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!

类别

在 Help Center 和 File Exchange 中查找有关 Contour Plots 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by