Adding a second y-axis with a different scale when fitting a curve.
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi everyone,
I am trying to add a second y-axis with a different scale using the attached code:
function ODE4522Aug2018 %-------------------------------------------------------------------------- function C=IntegratedModel(k,t)
c0 = [0.01721262399 0.00000008759 0.679173343 0.00000487921 0.00000487921 49.17075376105 0.00000000000];
options = odeset('RelTol',1e-3,'AbsTol',1e-3); [T,Cv]=ode45(@ODEloop,t,c0,options);
%disp('Table: Concentration of species '), disp (' c(1) c(2) c(3) (c4) (c5) (c6) c(7) k(1) k(2)') %disp([c(1)', c(2)', c(3)',c(4)',c(5)',c(6)', c(7)',k(1), k(2)'])
%fprintf('screen printing using fprintf\n') %fprintf('c(1)\t c(2)\t c(3)\t c(4)\t c(5)\t c(6)\t c(7)\t\n') %fprintf('%2.2f\t %5.5f\t %5.5f\t %5.5f\t %5.5f\t %5.5f\t %5.5f\t\n',[c(1);c(2);c(3);c(4);c(5);c(6);c(7)])
function dC=ODEloop(t,c)
%%PARAMETERS
A = 1.5e-6; B = 1.66667e-5; CC = 6.51332e-2; D = 0; E = 8.314; F = 323.15; G = 149; H = 4.14e-6; I = 1.39e-9; J = 2.89e-9; K = 8.4e-4; %k(2)= 9.598e-4; M = 5.15e+3; N = 3.53e-9; O = 1.07e-7; P = 10; Q = 8.825e-3; R = 12.54; S = 100.0869; %k(1)= 0.84; TT = 2703; U = 1.7e-3; V = 6.55e-8; W = 6.24; X =5.68e-5; Y = 258.30; Z = 2540; AA = 0.00000933254;
dcdt=zeros(7,1);
dcdt(1) = (1/A) * (B * CC - B * c(1)) - ((c(1) * E * F - G * (((c(3) * AA.^2) / (AA.^2 + W * AA + W * X))))/((1/H) + (G/((1 + (I* c(5))/(J* c(3))) * K))));
dcdt(2) = (1/A) * (B * D - B * c(2)) - (k(2) * (1 + (I* c(5))/(N * c(4)))) * (c(2)* E * F/M - (((c(3) * AA.^2) / (AA.^2 + U * AA + U * V))));
dcdt(3) = ((c(1) * E * F - G * ((c(3) * AA.^2) / (AA.^2 + W * AA + W * X)))/((1/H) + (G/((1 + (I* c(5))/(J* c(3)))*K)) - (0.162 * exp(-5153 / F) * ((( c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA + W * X) )) / O) - 1) * (P / (( c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA + W * X))) / O)))));
dcdt(4) = (k(2) * (1 + (I* c(5)) / (N * c(4))) * (((c(2)* E * F) / M)) - ((c(3) * AA.^2) / (AA.^2 + U * AA + U * V))) - (- Q * R * S * c(6) * c(7) * (1 - (k(1) * c(7)) / (1 + k(1) * c(7))));
dcdt(5) = (- Q * R * S * c(6) * c(7) *(1 -(k(1) * c(7)) / (1 + k(1) * c(7)))) - (0.162 * exp(- 5153/F) * (((c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA + W * X) )) / O) - 1) * (P / ((c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA + W * X) )) / O)));
dcdt(6) = - c(6) * (- Q * R * S * c(6) * c(7) * (1 - (k(1) * c(7)) / ( 1 + k(1) * c(7)))) * S / TT;
dcdt(7) = (- Q * R * S * c(6) * c(7) * (1 - (k(1) * c(7))/(1 + k(1) * c(7))))* (Y/ Z);
dC=dcdt; end C=Cv; end
t=[1200 2400 10200 ];
c=[ 0.01721262399 0.00000008759 0.679173343 0.00000487921 0.00000487921 49.17075376105 0.00000000000 0.01700907268 0.00000010281 1.405349320 0.00000511161 0.00000511161 49.60948155721 0.00000000000 0.01741617529 0.00000036495 10.714612593 0.00000535393 0.00000535393 30.91118867616 9.33482826985 ];
k0 = [0.3 9.598e-4];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@IntegratedModel,k0,t,c);
tv = linspace(min(t), max(t)); Cfit = IntegratedModel(k, tv);
figure(1) %plot(t, c, 'v') %hold on %hlp = plot(tv, Cfit); %hold off %xlabel('Time (sec)') %ylabel('Concentration (mol/m^3)') %legend(hlp, 'SO_{2,Headspace}', 'CO_{2,Headspace}','S_{total}','C_{total}','Ca^{2+}_{total}','CaCO_{3}', 'Ca^{2+}','Location','E')
yyaxis left plot(t, c(1),'rd',t, c(2),'cv',t, c(4),'y>',t, c(5),'m^') hold on hlp(1) = plot(tv, Cfit); hold off ylabel('Conc (mol/m^{3}') xlabel ('Time (sec)') set(gca,'YColor','k');
yyaxis right plot(t, c(3),'ks',t, c(6),'bo',c(7),'gp') hold on hlp(2) = plot(tv, Cfit); hold off legend(hlp(1),hlp(2), 'SO_{2,Headspace}', 'CO_{2,Headspace}','C_{total}','Ca^{2+}_{total}','S_{total}','CaCO_{3}', 'Ca^{2+}','location','Northeast') ylabel('Conc (mol/m^{3}') xlabel ('Time (sec)') set(gca,'YColor','k');
end
However, I get an error message that say:
"Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.
Error in ODE4522Aug2018 (line 112) hlp(1) = plot(tv, Cfit);"
What can I do?
2 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Gaussian Process Regression finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!