Error using lsqcurvefit (line 262) Function value and YDATA sizes are not equal.

Hi all,

I get the error message:

Error using lsqcurvefit (line 262)
Function value and YDATA sizes are not equal.

Error in IgorWaterODE15s (line 122)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

When I run this code:

function IgorWaterODE15s

    function C =kinetics(theta,t)
        
    % Initial conditions
c0 = [1; 1; 1; 1];

% Time
tspan = [0 30720];

% A constant, singular mass matrix
M = [1 0 0 0
   0 1 0 0
   0 0 0 0
   0 0 0 0];
% Options
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-6 1e-6 1e-6],'Vectorized','on');

% Solving
[t,Cv] = ode15s(@DifEq,tspan,c0,options);
%       
            
            function dC = DifEq(t,c)

            theta(1) = 5.97E-04;
            theta(2) =1.95E-03;
            theta(3) =8.314;
            theta(4) =323.15;
            theta(5) =5.3E-8;
            theta(6) =143.40;
            theta(7) =6.24;
            theta(8) =5.68E-05;
            theta(9) =2.04E-09;
            theta(10) =2.85E-09;
            theta(11) =1.70E-09;
            theta(12) =5.00E-06;
            theta(13) =4.72E-03;

            dC = [(theta(1)./theta(2)).*((1930.65./1000000)- (c(1,:) .*theta(3) .*theta(4) ./ 875000 )) - (c(1,:).* theta(3).* theta(4) - theta(6).* ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) - ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2) - ((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))).* theta(13)))
            (c(1,:).* theta(3).* theta(4) - theta(6).* ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) - ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2) - ((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))).* theta(13))) 
             - c(3,:) + (theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) + 2.* (theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2)+ (theta(5)./c(3,:))
            - c(4,:) + ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))) + 2.*((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))) + (theta(5)./c(4,:))];
            end
     
    C=Cv
    end



t=[0
960
1920
2880
3840
4800
5760
6720
7680
8640
9600
10560
11520
12480
13440
14400
15360
16320
17280
18240
19200
20160
21120
22080
23040
24000
24960
25920
26880
27840
28800
29760
30720];


c=  [3.16E-02	0.33	1.62E-04	1.35E-04
    3.56E-02	0.64	2.69E-04	2.24E-04
    3.83E-02	0.93	4.46E-04	3.72E-04
    4.12E-02	1.18	2.51E-03	2.09E-03
    4.34E-02	1.40	0.43        0.35
    4.58E-02	1.60	1.20        1.00
    4.74E-02	1.78	1.70        1.41
    4.92E-02	1.94	1.99        1.66
    5.08E-02	2.08	2.81        2.34
    5.21E-02	2.20	2.95        2.45
    5.34E-02	2.31	3.23        2.69
    5.49E-02	2.41	3.31        2.75
    5.56E-02	2.49	3.38        2.82
    5.63E-02	2.57	3.16        2.63
    5.72E-02	2.63	3.08        2.57
    5.80E-02	2.69	3.38        2.82
    5.88E-02	2.74	3.23        2.69
    5.92E-02	2.78	3.23        2.69
    5.98E-02	2.82	3.62        3.02
    6.03E-02	2.85	3.71        3.09
    6.06E-02	2.87	3.71        3.09
    6.10E-02	2.89	3.54        2.95
    6.14E-02	2.91	3.54        2.95
    6.15E-02	2.93	3.71        3.09
    6.18E-02	2.94	3.79        3.16
    6.20E-02	2.95	3.71        3.09
    6.22E-02	2.96	3.88        3.24
    6.24E-02	2.96	3.88        3.24
    6.25E-02	2.97	3.88        3.24
    6.26E-02	2.97	4.07        3.39
    6.27E-02	2.97	4.07        3.39
    6.28E-02	2.98	4.16        3.47
    6.29E-02	2.98	4.16        3.47]; 


theta0=[1;1;1;1;1;1;1;1;1;1;1;1;1];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);


fprintf(1,'\tRate Constants:\n')
 for k1 = 1:length(theta)
     fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
 end

tv = linspace(min(t), max(t),33);
tv = tv';


Cv = kinetics(theta, tv)

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cv);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

end

Where am I going wrong?

Answer 1

[Torsten]

In "kinetics", use

tspan = t;

instead of

tspan = [0 30720];

[Dursman Mchabe]

Thanks a lot Torsten. It works.

Answer 2

[Walter Roberson]

You have

function C =kinetics(theta,t)

The input values in the first parameter, theta, will be the proposed parameter values. The second parameter, t, will be the "x" values at which the function needs to calculate values given the input proposed parameter values and the x values. There must be one output value for each "x" value, so the output must C must have the same number of elements as t.

Inside the function you have

tspan = [0 30720];
[t,Cv] = ode15s(@DifEq,tspan,c0,options);
C=Cv;

You have ignored the input t values (and so ignored how many there are), and you overwrite the variable t with the times at which projected values for DifEq were calculated.

The number of output t values will be the same as the number of output rows in Cv -- each Cv row gives response at the corresponding (output) t value.

The ode* routines are variable step solvers and adjust the number of locations they actually evaluate at in order to meet integration targets. However, they do output exactly once per time internally.Instead they work in one of two modes:

1. if tspan was passed in with three or more elements, then the output times are the same as the times in tspan, but the times internally evaluated at could be quite different, either more (if the function is oscillating badly) or fewer (a very smooth function might not need to be evaluated at a particular time in order to project the value within target accuracy.)
2. if tspan was passed in with exactly two elements, as in your case, then the output times are determined internally, and not all times evaluated at are necessarily reported as they might not turn out to add much information about the behavior of the system; furthermore, the algorithm feels free to report at times that it did not specifically evaluate it (the times reported at are synthesized to keep the errors within bounds according to the projections.)

If what you are doing is effectively boundary value problems, then you should (A) probably be using a bvp* routine instead of what you are doing now; and (B) you should probably be passing in the input t as the tspan array.

[Dursman Mchabe]

Thanks a lot Walter, when I use t, instead of tspan, the code runs.