% Plotting code for the PCM_CODE.m. % This will plot three different times in separate windows for times that are % specified. % % Cameron Smith % Date Created: 05/02/16 % Last Modified: 25/02/16 close all clear % Firstly specify the possible intial conditions % Set this to one if the initial mass is in the left half (PDE half) of the domain INITIAL_LEFT = 1; % Set this to one if the initial mass is in the right half (compartment half) of the domain INITIAL_RIGHT = 0; % Set this to one if you want a uniform initial condition INITIAL_UNIFORM = 0; % Set this to one if you want all of the mass in a compartments width starting at the lower domain end INITIAL_COMP = 0; % Sum these intial logicals SUM = INITIAL_LEFT+INITIAL_RIGHT+INITIAL_UNIFORM+INITIAL_COMP; % Check that we have only one initial condition if (SUM > 1) error('Please choose only one initial condition, you have chosen %.0f',SUM) elseif (SUM < 1) error('Please choose an initial condition') end fprintf('Loading Workspace...') % Load the correct workspace if INITIAL_LEFT == 1 load('WORKSPACE_PCM_INITIAL_LEFT') elseif INITIAL_RIGHT == 1 load('WORKSPACE_PCM_INITIAL_RIGHT') elseif INITIAL_UNIFORM == 1 load('WORKSPACE_PCM_INITIAL_UNIFORM') elseif INITIAL_COMP == 1 load('WORKSPACE_PCM_INITIAL_COMP') end fprintf('DONE\n') % Create the plotting ranges x_PLOT_PDE = x_MIN:h_PDE:x_MAX; x_PLOT_COM = x_MIN+h_COM/2:h_COM:x_MAX-h_COM/2; % Create the recording vector t_REC_VEC = t_INITIAL:t_REC_STEP:t_FINAL; % Specify a vector of times to plot at t_PLOT_TIMES = [t_INITIAL,t_FINAL/4,t_FINAL]; % Find the corresponding indices for the PDE and compartment regimes t_PLOT_PDE_INDEX = ceil(t_PLOT_TIMES./t_FINAL.*(t_PDE_LEN-1))+1; t_PLOT_COM_INDEX = ceil(t_PLOT_TIMES./t_FINAL.*(t_REC_LEN-1))+1; % Plotting for i = 1:3 % Start a new figure figure % Hold the plots hold on % Extract indices IND_PDE = t_PLOT_PDE_INDEX(i); IND_COM = t_PLOT_COM_INDEX(i); % Plot the compartment solution in blue bars bar(x_VEC_COM_MID,SOL_COM_MEAN(IND_COM,:)/h_COM) % Plot the PDE solution in green plot(x_VEC_PDE,SOL_PDE_MEAN(IND_COM,:),'g','LineWidth',2) % Plot the interface as a red line plot([I,I],[0,max(max(SOL_PDE_MEAN(IND_COM,:)),max(SOL_COM_MEAN(IND_COM,:)/h_COM))+10],'r','LineWidth',3) % Plot the numerical solution plot(x_PLOT_PDE,SOL_PDE_EX(IND_PDE,:),'k--','LineWidth',2) % Add some labels titlestr = sprintf('Hybrid method for time t = %.2f',t_PLOT_TIMES(i)); title(titlestr) xlabel('x') end