Filter löschen
Filter löschen

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

1 Ansicht (letzte 30 Tage)
Dursman Mchabe
Dursman Mchabe am 22 Aug. 2018
Kommentiert: Dursman Mchabe am 22 Aug. 2018
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
jonas
jonas am 22 Aug. 2018
You can tell us the size and class of tv and Cfit
Dursman Mchabe
Dursman Mchabe am 22 Aug. 2018
tv size = t size (3 x 1) and Cv size = c size (3 x 7)

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Kategorien

Mehr zu Programming 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!

Translated by