Using yyaxis with lsqcurvefit

Hi everyone, I'm trying to use yyaxis with lsqcurve fit on this 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 ylim([0,0.1]) 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 ylim([0,50]) 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

Yet I get the error message:

"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 try/do?

[dpb]

Remove subscripting on hlp, set breakpoint and see what is returned. It appears plot must be returning a vector of handles rather than just one...

[Dursman Mchabe]

Thanks a lot for the response. When I remove subscripting on hlp I get the error message :

" Error using plot Data must be a single matrix Y or a list of pairs X,Y.

Error in ODE4522Aug2018 (line 120) plot(t, c(3),'ks',t, c(6),'bo',c(7),'gp')"

The plot comes out as the figure on the attachment.

[Walter Roberson]

Your t is a vector but c(3) is a scalar

[Dursman Mchabe]

Thanks a lot for the comment Walter. I have tried to convert c to vector, using vec=ind2vec(c).

And it seems I'm not licensed, I got the error message: "To use 'ind2vec', the following product must be licensed, installed, and enabled: Neural Network Toolbox

Error in ODE4522Aug2018 (line 92) vec=ind2vec(c);"

Is there any way around it?

[dpb]

What are you actually trying to plot? Why are you calling plot with the vector t for X but only a single value for the Y, c(3)?

If you expand c(3) to size/length of t, you'll just get a horizontal line of markers; is that the intent? If so, then the syntax would be

c(3)*ones(size(t))

to generate a vector the size of t with value c(3). Of course, you've got the same issue with the other values in the plot statement so you've got to fix 'em all (or change what it is you're trying to do, we can't know that).

nT=size(t);
plot(t, c(3)*ones(nT),'ks',t, c(6)*ones(nT),'bo',c(7)*ones(nT),'gp')

[Dursman Mchabe]

Thanks a lot for pointing me in the right direction. I need to generate a vector size of t for c(1), c(2), c(3), c(4), c(5), c(6) and c(7). When I try:

nT=size(t);
plot(t, c(3)*ones(nT),'ks',t, c(6)*ones(nT),'bo',c(7)*ones(nT),'gp'),

The entire plot becomes blank. I don't know how to do it.

[dpb]

Well, it "works" for whatever it is worth if there are data...

clear c
c(3)=pi; c(6)=8;
figure
plot(t,c(3)*ones(nT),'ks',t,c(6)*ones(nT),'bo')
xlim([0 11]),ylim([0 10])

gives

which illustrates the code snippet itself works. Have you verified that all c() are not NaN? plot will just ignore NaN silently.

[Walter Roberson]

Perhaps you want to use refline() or vertline()

[Dursman Mchabe]

Thanks a lot for the help. I truly appreciate. It works.

[dpb]

So, what was the end result of what was the problem? (inquiring minds and all that... :) )

[Dursman Mchabe]

Hello dpb, sorry for the late response. I was away for a weekend. The obtained graph is as on the attachment. I think the challenge is to express Cfit in terms of Cfit(1), Cfit(2)... Cfit(7).