Plot coming up blank

3 vues (au cours des 30 derniers jours)
Dursman Mchabe
Dursman Mchabe le 20 Août 2018
Modifié(e) : Walter Roberson le 21 Août 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.

Connectez-vous pour répondre à cette question.

Réponses (1)

Walter Roberson
Walter Roberson le 20 Août 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 commentaires
Dursman Mchabe
Dursman Mchabe le 21 Août 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 le 21 Août 2018
Modifié(e) : Walter Roberson le 21 Août 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,:))];

Connectez-vous pour commenter.

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by