How to stop matlab using previously calculated values in a subsequent calculation

Hello,

I am carrying out calculations using three equations for differing time steps. However when plot the results for the third time step it is not giving me the correct graph. I know this as I have commented the first and second time step calculations and it is giving me the correct plot. I assume it is because the programme is using previously calculated values but I cannot figure out which ones and how to stop it. Any help would be greatly appreciated. Thanks!

%% Clear workspace

clear all;

clc;

%% Input data

u= 5.5; % Mean wind velocity (m/s)

L= 10000; % Length of air quality concentration model (m)

dx= [500,900,1000,1100,1500];% step in x direction

dt= dx./u; % time step in seconds

Nt= ((10*3600)./dt)+1; % number of time steps

Nj= 10; % number of spacial steps

c0= 0; % No emission concentration (mg/m^3)

nstop=9.5*3600./dt;

dthours=dt./3600;

tmin=0;

tmax=10; % hours

T= tmin:dthours:tmax; %time period

clb=93.6; % Emission of SO2 (mg/m^3)

%% Allocate memory for matricies of differing time steps

%LAX

c1=zeros(Nt(1),Nj+1);

c2=zeros(Nt(2),Nj+1);

c3=zeros(Nt(3),Nj+1);

c4=zeros(Nt(4),Nj+1);

c5=zeros(Nt(5),Nj+1);

%FTCS

cf1=zeros(Nt(1),Nj+1);

cf2=zeros(Nt(2),Nj+1);

cf3=zeros(Nt(3),Nj+1);

cf4=zeros(Nt(4),Nj+1);

cf5=zeros(Nt(5),Nj+1);

%Upwind

cu1=zeros(Nt(1),Nj+1);

cu2=zeros(Nt(2),Nj+1);

cu3=zeros(Nt(3),Nj+1);

cu4=zeros(Nt(4),Nj+1);

cu5=zeros(Nt(5),Nj+1);

%% Initial condition

c1(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

c2(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

c3(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

c4(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

c5(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cf1(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cf2(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cf3(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cf4(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cf5(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cu1(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cu2(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cu3(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cu4(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

cu5(1,:)=0.0011; % Initial concentration of SO2 (mg/m^3)

%% Main calculation

% dx and dt (1)

for n=1:Nt(1)-1

if n < nstop(1)

c1(n,1)=clb;

cf1(n,1)=clb;

cu1(n,1)=clb;

else

clb=c0;

end

for j=2:Nj

c1(n+1,j) = c1(n,j-1); %LAX Scheme

cf1(n+1,j) = cf1(n,j)-0.5 * cf1(n,j+1)+0.5 * cf1(n,j-1); %FTCS Scheme

cu1(n+1,j) = cu1(n,j-1); %Upwind Scheme

end

c1(n+1,Nj+1)=c1(n+1,Nj);

cf1(n+1,Nj+1)=cf1(n+1,Nj);

cu1(n+1,Nj+1)=cu1(n+1,Nj);

end

clear n

% dx and dt (2)

for n=1:Nt(2)-1

if n < nstop(2)

c2(n,1)=clb;

cf2(n,1)=clb;

cu2(n,1)=clb;

else

clb=c0;

end

for j=2:Nj

c2(n+1,j) = c2(n,j-1); %LAX Scheme

cf2(n+1,j) = cf2(n,j)-0.5 * cf2(n,j+1)+0.5 * cf2(n,j-1); %FTCS Scheme

cu2(n+1,j) = cu2(n,j-1); %Upwind Scheme

end

c2(n+1,Nj+1)=c2(n+1,Nj);

cf2(n+1,Nj+1)=cf2(n+1,Nj);

cu2(n+1,Nj+1)=cu2(n+1,Nj);

end

clear n

% dx and dt (3)

for n=1:Nt(3)-1

if n < nstop(3)

c3(n,1)=clb;

cf3(n,1)=clb;

cu3(n,1)=clb;

else

clb=c0;

end

for j=2:Nj

c3(n+1,j) = c3(n,j-1); %LAX Scheme

cf3(n+1,j) = cf3(n,j)-0.5 * cf3(n,j+1)+0.5 * cf3(n,j-1); %FTCS Scheme

cu3(n+1,j) = cu3(n,j-1); %Upwind Scheme

end

c3(n+1,Nj+1)=c3(n+1,Nj);

cf3(n+1,Nj+1)=cf3(n+1,Nj);

cu3(n+1,Nj+1)=cu3(n+1,Nj);

end

%% Plot for timestep itteration 3

elax3=c3(end-11:end-1,:);

slax3=c3(end-(198):end-189,:);

eFTCS3=cf3(end-11:end-1,:);

sFTCS3=cf3(end-198:end-189,:);

eUpwind3=cu3(end-11:end-1,:);

sUpwind3=cu3(end-198:end-189,:);

disp('Startup Concentraion Profile (LAX SCHEME):')

disp(slax3)

disp('Shutdown Concentraion Profile (LAX SCHEME):')

disp(elax3)

disp('Startup Concentraion Profile (FTCS SCHEME):')

disp(sFTCS3)

disp('Shutdown Concentraion Profile (FTCS SCHEME):')

disp(eFTCS3)

disp('Startup Concentraion Profile (UPWIND SCHEME):')

disp(sUpwind3)

disp('Shutdown Concentraion Profile (UPWIND SCHEME):')

disp(eUpwind3)

figure(2)

hold on

plot(slax3(:,3),'blue')

hold on

plot(sFTCS3(:,3),'red')

hold on

plot(sUpwind3(:,3),'green','linestyle','--')

xlabel('Time (Hours)')

ylabel('Concentration (mg/m^3)' )

Answer 1

[Maximilian]

Problem has been solved thank you