Adding a second y-axis with a different scale when fitting a curve.

1 次查看(过去 30 天)
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?

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by