Adding a second y-axis with a different scale when fitting a curve.
1 view (last 30 days)
Show older comments
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 Comments
Answers (0)
See Also
Categories
Find more on Acquisition Using Kinect for Windows Hardware in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!