Plot coming up blank

Question

Dursman Mchabe on 20 Aug 2018

0
Link

Direct link to this question

https://ch.mathworks.com/matlabcentral/answers/415579-plot-coming-up-blank

Edited: Walter Roberson on 21 Aug 2018

ODE15s20Aug2018.m

Open in MATLAB Online

Hi everyone,

I am trying to plot

function  ODE15s20Aug2018
%%VARIABLES
% y(1)
% y(2)
% y(3) 
% y(4)
% y(5)
% y(6) 
% y(7) 
% y(8) 
%%EQUATIONS 
%(y(1)' = ((1/A) * (B * C– B * y(1))) – (((y(1) * E * F) – (G * ((y(3) * y(8).^2)/(y(8).^2 + H * y(8) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5)) / (LL * y(3)))) * M))))
%(y(2)' = ((1/A) * (B * N– B * y(2))) – ((O * (1 + ((K * y(5))/ (P * y(4)))))* (((y(2) * E * F) / (S)) – (((y(4) * y(8).^2) / (y(8).^2 + Q * y(8) + Q * R)))))
%(y(3)) /dt = ( ((y(1) * E * F) – (G * ((y(3) * y(8).^2)/(y(8).^2 + H * y(8) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5)) / (LL * y(3)))) * M)))) – (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T))))
%(y(4)) /dt = ((O* (1 + ((K * y(5))/ (P * y(4)))))* (((y(2) * E * F) / (S)) – (((y(4) * y(8).^2) / (y(8).^2 + Q * y(8) + Q * R))))) + (y(6) * (V* (1 – ((y(5) * ((y(4) *  Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W))))
%(y(5)' = (y(6) * (V* (1 – ((y(5) * ((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W)))) – (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T))))
%(y(6)' = -y(6) * (V* (1 – ((y(5) * ((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W))) * (X / Y)
%(y(7)' = (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T)))) * (Z/AA)
%y(8) + 2 * y(5) – ((y(3) * H * y(8))/(y(8).^2 + H * y(8) + H * I)) – 2 * (((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) – (((y(4) * Q * y(8))/ (y(8).^2 + Q * y(8) + Q * R))) – 2 * (((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) – (AB/y(8)) = 0
%%SINGULAR MASS MATRIX
M = diag([1 1 1 1 1 1 1 0 ]);
%%INITIAL VALUES 
y0 = zeros(8,1);
y0(1)= 1e-9; 
y0(2)= 1e-9; 
y0(3)= 1e-9; 
y0(4)= 1e-9; 
y0(5)= 1e-9; 
y0(6)= 4.99e-9; 
y0(7)= 1e-9; 
y0(8)= 7.413e-6; 
tspan = linspace(0,183000,30);
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-4,'Vectorized','off');
[t,y] = ode15s(@f,tspan,y0,options);
plot(t,y)
% --------------------------------------------------------------------------
function out = f(t,y)
% PARAMETERS
A = 1.5e-6;      
B = 1.67e-5;    
C = 6.51e-2;    
E = 8.314;      
F = 323.15;     
G = 149;        
H = 6.24;       
I = 5.68e-5;    
J = 4.14e-6;    
K = 1.39e-9;    
LL = 2.89e-9;   
M = 8.4e-4;     
N = 0;          
O = 9.598e-4;   
P = 3.53e-9;    
Q = 1.7e-3;    
R = 6.55e-8;    
S = 5.15e3;     
T = 1.07e-7;    
U = 10;         
V = 8.42e-4;     
W = 3.2e-1;     
X = 100.086;    
Y = 2703;       
Z = 258.3;       
AA = 2540;      
AB = 5.3e-8;    
   out =[ ((1/A) * (B * C - B * y(1,:))) - (((y(1,:) * E * F) - (G * ((y(3,:) * y(8,:).^2)/(y(8,:).^2 + H * y(8,:) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5,:)) / (LL * y(3,:)))) * M))))
      ((1/A) * (B * N - B * y(2,:))) - ((O * (1 + ((K * y(5,:))/ (P * y(4,:)))))* (((y(2,:) * E * F) / (S)) - (((y(4,:) * y(8,:).^2) / (y(8,:).^2 + Q * y(8,:) + Q * R)))))
      ( ((y(1,:) * E * F) - (G * ((y(3,:) * y(8,:).^2)/(y(8,:).^2 + H * y(8,:) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5,:)) / (LL * y(3,:)))) * M)))) - (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T))))
      ((O* (1 + ((K * y(5,:))/ (P * y(4,:)))))* (((y(2,:) * E * F) / (S)) - (((y(4,:) * y(8,:).^2) / (y(8,:).^2 + Q * y(8,:) + Q * R))))) + (y(6,:) * (V* (1 - ((y(5,:) * ((y(4,:) *  Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W))))
      (y(6,:) * (V* (1 -  ((y(5,:) * ((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W)))) - (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T))))
      -y(6,:) * (V* (1 - ((y(5,:) * ((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W))) * (X / Y)
      (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T)))) * (Z/AA)
      y(8,:) + 2 * y(5,:) - ((y(3,:) * H * y(8,:))/(y(8,:).^2 + H * y(8,:) + H * I)) - 2 * (((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) - (((y(4,:) * Q * y(8,:))/ (y(8,:).^2 + Q * y(8,:) + Q * R))) - 2 * (((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) - (AB/y(8,:))];

However, my plot comes out blank. What can I do?

See attached for details.

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Sign in to answer this question.

Answer 1

Walter Roberson on 20 Aug 2018

0
Link

Direct link to this answer

https://ch.mathworks.com/matlabcentral/answers/415579-plot-coming-up-blank#answer_333480

Open in MATLAB Online

>> ODE15s20Aug2018
Warning: Failure at t=5.450176e-06.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.355253e-20) at time t. 
> In ode15s (line 668)
  In ODE15s20Aug2018 (line 48)

Your function contains a singularity that is encountered in the very first time step after time 0, so your t is coming out as a single value, and your y is coming out as a single row vector. When you ask to plot() a single point, then because by default it only plots lines, you do not see anything in your plot. If you were to add a marker style in your plot, you would have gotten a plot with a single point drawn.

If you have the symbolic toolbox, I recommend that you look at the documentation for odeFunction() and follow the first example there to see the workflow of converting a symbolic system of equations for use with the numeric ode* routines.

2 Comments
Show NoneHide None

Dursman Mchabe on 21 Aug 2018

Thanks Walter, I will check the documentation. My end-goal is to compare y(1) to y(8) with the experimental results that were collected over a period of 183 000 seconds.

Dursman Mchabe on 21 Aug 2018

Edited: Walter Roberson on 21 Aug 2018

Open in MATLAB Online

function  ODE15s19Aug2018
%%VARIABLES
% y(1) = CSO2_Headspace
% y(2) = CCO2_Headspace
% y(3) = S_total
% y(4) = C_total
% y(5) = CCa2
% y(6) = CCaCO3
% y(7) = CCaSO3
% y(8) = CH
%%EQUATIONS 
%(y(1)' = ((1/A) * (B * C– B * y(1))) – (((y(1) * E * F) – (G * ((y(3) * y(8).^2)/(y(8).^2 + H * y(8) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5)) / (LL * y(3)))) * M))))
%(y(2)' = ((1/A) * (B * N– B * y(2))) – ((O * (1 + ((K * y(5))/ (P * y(4)))))* (((y(2) * E * F) / (S)) – (((y(4) * y(8).^2) / (y(8).^2 + Q * y(8) + Q * R)))))
%(y(3)) /dt = ( ((y(1) * E * F) – (G * ((y(3) * y(8).^2)/(y(8).^2 + H * y(8) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5)) / (LL * y(3)))) * M)))) – (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T))))
%(y(4)) /dt = ((O* (1 + ((K * y(5))/ (P * y(4)))))* (((y(2) * E * F) / (S)) – (((y(4) * y(8).^2) / (y(8).^2 + Q * y(8) + Q * R))))) + (y(6) * (V* (1 – ((y(5) * ((y(4) *  Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W))))
%(y(5)' = (y(6) * (V* (1 – ((y(5) * ((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W)))) – (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T))))
%(y(6)' = -y(6) * (V* (1 – ((y(5) * ((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) / W))) * (X / Y)
%(y(7)' = (0.162 * (exp(-5153/F)) * (((y(5) * ((y(3) * H * I-) / (y(8).^2 + H * y(8) + H * I))) / T) - 1).^3 * ( (U) / (( y(5) * ((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) /(T)))) * (Z/AA)
%y(8) + 2 * y(5) – ((y(3) * H * y(8))/(y(8).^2 + H * y(8) + H * I)) – 2 * (((y(3) * H * I) / (y(8).^2 + H * y(8) + H * I))) – (((y(4) * Q * y(8))/ (y(8).^2 + Q * y(8) + Q * R))) – 2 * (((y(4) * Q * R) / (y(8).^2 + Q * y(8) + Q * R))) – (AB/y(8)) = 0
%%SINGULAR MASS MATRIX
M = diag([1 1 1 1 1 1 1 0 ]);
%%INITIAL VALUES 
y0 = zeros(8,1);
y0(1)= 1e-9; 
y0(2)= 1e-9; 
y0(3)= 1e-9; 
y0(4)= 1e-9; 
y0(5)= 1e-9; 
y0(6)= 4.99e-9; 
y0(7)= 1e-9; 
y0(8)= 7.413e-6; 
tspan = linspace(0,183000,30);
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-4,'Vectorized','off');
[t,y] = ode15s(@f,tspan,y0,options);
plot(t,y)
%y(8,:) = -log(y(8,:));
%%EXPERIMENTAL RESULTS
%Time = [600;1200;1800;2400;3000;10200;17400;24600;31800;39000;46200;53400;60600;67800;75000;82200;89400;...
%96600;103800;111000;118200;125400;132600;139800;147000;154200;161400;168600;175800;183000;];
%SO2Headspace = [0.017599209;0.017212624;0.017111011;0.017009073;0.016805847;0.017416175;0.016500683;...
%0.016704234;0.016704234;0.017009073;0.018005986;0.021504463;0.026813406;0.030718986;0.033334539;0.036841483;...
%0.041438486;0.044855868;0.048659509;0.051954761;0.054761815;0.057568868;0.059765594;0.061779612;0.06336601;...
%0.064484728;0.064871312;0.065379702;0.065684866;0.065888417];
%S_total = [3.2937e1;6.7022e1;1.0226e2;1.3864e2;5.6500e2;1.0023e3;1.4582e3;1.9271e3;2.4127e3;...
%2.9215e3;3.4154e3;3.8883e3;4.3330e3;4.7362e3;5.1279e3;5.4995e3;5.8130e3;6.1206e3;...
%6.4217e3;6.6661e3;6.9004e3;7.1184e3;7.3118e3;7.4887e3;7.6451e3;7.7850e3;7.7971e3;...
%7.8030e3;7.8053e3;7.8052e3];
%Ca2 = [0.004879211;0.004879211;0.005111607;0.005111607;0.005353935;0.005353935;0.005363117;0.005363117;...
%0.005073082;0.005073082;0.004933355;0.005296472;0.005296472;0.005296472;0.009456959;0.009456959;0.011300764;...
%0.011457084;0.01099768;0.01099768;0.010393233;0.01038722;0.01038722;0.010256949;0.010364689;0.01055387;...
%0.011318728;0.011318728;0.011318728;0.010264908;];
%H = [7.4131E-06;9.33254E-06;1.20226E-05;1.04713E-05;1.7378E-05;3.71535E-05;3.54813E-05;0.0001;...
%0.040738028;0.245470892;0.616595002;1.318256739;2.290867653;2.344228815;2.398832919;1.819700859;...
%1.380384265;2.089296131;1.819700859;1.584893192;2.290867653;1.621810097;1.122018454;0.891250938;...
%0.851138038;0.758577575;0.87096359;1.122018454;1;0.851138038;];
%CaCO3 = [48.62382123;49.17075376;49.43940249;49.60948156;37.76118391;30.91118868;26.72559345;19.77300722;...
%11.40329066;9.577469036;5.65344167;3.315502675;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;];
%CaSO3 = [0.0000;0.0000;0.0000;0.0000;6.3007;9.3348;11.2637;14.2916;17.8975;19.0564;20.8685;22.0735;...
%23.7276;24.0230;24.3478;24.7041;24.9783;25.3534;25.8385;26.2147;26.6544;27.1365;27.6080;28.0962;28.5717;...
%29.0324;29.0324;29.0324;29.0324;29.0324;];
%pH = [8.13;8.03;7.92;7.98;7.76;7.43;7.45;7;4.39;3.61;3.21;2.88;2.64;2.63;2.62;2.74;...
%2.86;2.68;2.74;2.8;2.64;2.79;2.95;3.05;3.07;3.12;3.06;2.95;3;3.07;];
%figure(1)
%yyaxis left          
%plot(Time,CaCO3,'kd',t,y(6,:),'k-.',Time,CaSO3,'m*',Time,pH,'r>',t,y(8,:),'r-.')
%plot(t,y(:,6),'k-.',t,y(:,8),'r-.')
%ylabel('Conc (mol/m^{3}')
%xlabel ('Time (sec)')
%set(gca,'YColor','k');
%yyaxis right
%plot(Time,SO2Headspace,'bo',t,y(1,:),'b-')
%plot(t,y(:,1),'b-')
%legend('Exp CaCO_{3}','Mod CaCO_{3}','Exp CaSO_{3}.^{1}/_{2}H_{2}O','Exp pH','Mod pH',...
%'Exp SO_{2,Headspace}','Mod SO_{2,Headspace}','location','Northeast')
%ylabel('Conc (mol/m^{3}')
%xlabel ('Time (sec)')
%set(gca,'YColor','k');
% --------------------------------------------------------------------------
function out = f(t,y)
% PARAMETERS
A = 1.5e-6;     % V_Headspace (m3)     Volume of headspace 
B = 1.67e-5;    % F      (m3/s)   Flue gas flow rate
C = 6.51e-2;    % CSO2_in    (mol/m3)  Inlet flue gas SO2 concentration
E = 8.314;      % R      (J/molK)  Universal gas constant
F = 323.15;     % T       (K)    Reactor temperature
G = 149;        % HSO2    (Pam3/mol)Henry’s constant
H = 6.24;       % KSO2    (mol/m3)  SO2 dissociation equilibrium constant
I = 5.68e-5;    % KHSO3    (mol/m3)  HSO3- dissociation equilibrium constant
J = 4.14e-6;    % kga      (1/s)  Gas-side mass transfer-surface area SO2
K = 1.39e-9;    % DCa2    (m2/s)  Diffusivity of calcium ions
LL = 2.89e-9;   % DSO2    (m2/s)  Diffusivity of SO2 
M = 8.4e-4;     % kLa      (1/s)  Liquid-side mass transfer-surface area SO2
N = 0;          %CCO2_in    (mol/m3)  Inlet flue gas CO2 concentration
O = 9.598e-4;   % kLa_CO2    (1/s)  Liquid-side mass transfer-surface area CO2 
P = 3.53e-9;    % DCO2    (m2/s)  Diffusivity of CO2
Q = 1.7e-3;     % KCO2    (mol/m3)  CO2 dissociation equilibrium constant
R = 6.55e-8;    % KHCO3    (mol/m3)  HCO3- dissociation equilibrium constant
S = 5.15e3;     % HCO2    (Pam3/mol)CO2 Henry’s law constant
T = 1.07e-7;    % KSPCaSO3    (mol2/m6)  Solubility product of CaSO3 
U = 10;         % BETCaSO3    (m2/g)  BET specific surface area
V = 8.42e-4;     % kd_CaCO3    (kg/mol.s)CaCO3 dissolution constant 
W = 3.2e-1;     % KSP,CaCO3    (mol2/m6)  Solubility product of CaCO3
X = 100.086;    % MWCaCO3    (g/mol)  Molecular weight of CaCO3
Y = 2703;       % rhoCaCO3    (kg/m3)  Density of CaCO3
Z = 258.3;       % MWCaSO3    (g/mol)  Molecular weight of CaSO3.0.5H2O
AA = 2540;      % rhoCaSO3    (kg/m3)  Density of CaSO3.0.5H2O
AB = 5.3e-8;    % Kw      (mol2/m6) H2O dissociation equilibrium constant
   out =[ ((1/A) * (B * C - B * y(1,:))) - (((y(1,:) * E * F) - (G * ((y(3,:) * y(8,:).^2)/(y(8,:).^2 + H * y(8,:) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5,:)) / (LL * y(3,:)))) * M))))
      ((1/A) * (B * N - B * y(2,:))) - ((O * (1 + ((K * y(5,:))/ (P * y(4,:)))))* (((y(2,:) * E * F) / (S)) - (((y(4,:) * y(8,:).^2) / (y(8,:).^2 + Q * y(8,:) + Q * R)))))
      ( ((y(1,:) * E * F) - (G * ((y(3,:) * y(8,:).^2)/(y(8,:).^2 + H * y(8,:) + H * I))))/ ((1 / J) + (G / ((1+ ((K * y(5,:)) / (LL * y(3,:)))) * M)))) - (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T))))
      ((O* (1 + ((K * y(5,:))/ (P * y(4,:)))))* (((y(2,:) * E * F) / (S)) - (((y(4,:) * y(8,:).^2) / (y(8,:).^2 + Q * y(8,:) + Q * R))))) + (y(6,:) * (V* (1 - ((y(5,:) * ((y(4,:) *  Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W))))
      (y(6,:) * (V* (1 -  ((y(5,:) * ((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W)))) - (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T))))
      -y(6,:) * (V* (1 - ((y(5,:) * ((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) / W))) * (X / Y)
      (0.162 * (exp(-5153/F)) * (((y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) / T) - 1).^3 * ( (U) / (( y(5,:) * ((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) /(T)))) * (Z/AA)
      y(8,:) + 2 * y(5,:) - ((y(3,:) * H * y(8,:))/(y(8,:).^2 + H * y(8,:) + H * I)) - 2 * (((y(3,:) * H * I) / (y(8,:).^2 + H * y(8,:) + H * I))) - (((y(4,:) * Q * y(8,:))/ (y(8,:).^2 + Q * y(8,:) + Q * R))) - 2 * (((y(4,:) * Q * R) / (y(8,:).^2 + Q * y(8,:) + Q * R))) - (AB/y(8,:))];

Sign in to comment.

Plot coming up blank

0 Comments
Show -2 older commentsHide -2 older comments

Answers (1)

2 Comments
Show NoneHide None

See Also

Categories

Tags

Community Treasure Hunt

Plot coming up blank

0 Comments Show -2 older commentsHide -2 older comments

Answers (1)

2 Comments Show NoneHide None

See Also

Categories

Tags

Community Treasure Hunt

0 Comments
Show -2 older commentsHide -2 older comments

2 Comments
Show NoneHide None