Plot coming up blank

3 visualizzazioni (ultimi 30 giorni)
Dursman Mchabe
Dursman Mchabe il 20 Ago 2018
Modificato: Walter Roberson il 21 Ago 2018
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.

Accedi per rispondere a questa domanda.

Risposte (1)

Walter Roberson
Walter Roberson il 20 Ago 2018
>> 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 Commenti
Dursman Mchabe
Dursman Mchabe il 21 Ago 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
Dursman Mchabe il 21 Ago 2018
Modificato: Walter Roberson il 21 Ago 2018
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,:))];

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by