Fitting of two differential equations simultaneously using lsqcurvefit

Hello,

I am new to matlab and I have some questions regarding fitting of two curves simultaneously using lsqcurvefit. Currently, I have two sets of experimental datas (let's say A and B) and would like to fit into my model that has two differential equations (dA/dt, dB/dt) with two parameters (B(1) and B(2)) I have read the previous post and wrote the code as:

[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel, B0, t(:), [A B], [0;0], [Inf;Inf]);

with plmodel consist of two differential equations (in the form of "ode45()"), B0 = initial guess of the parameters (randomly generated), t(:) = input (in (1xn) matrix), [A B] = experimental data.

The problem I am facing is that it seems to fit one of the data only (almost perfectly), while the other data is totally off. It also show low resnorm (0.0266) for this fitting. The questions I would want to ask are:

1. Why does the matlab give good resnorm for the fitting eventhough one of the data does not fit very well?
2. Is there any way to fix this so that I can fit both data simultaneously?

Thank you!

[Torsten]

Code ?

[Solomon Thong]

Hi, this is the code. Thank you:

function curveFitting
%Data, t = time; Aconc = concentration of A; Bconc = concentration of B
t = [0 5 10 20 30 45 60];
Aconc = 8*[1;0.546547;0.339339;0.215465;0.175676;0.15991;0.153153];
Cconc = 0.993*[1;0.683173;0.556456;0.448287;0.382121;0.316049;0.257132];
 
 function S = plmodel(B,t) 
 x0 = [8;0.0993]; %%Initial concentration
 [T,Sv] = ode45(@DifEq,t,x0); %%Runge Kutta 4th order
     function dS = DifEq(t,x)  
     xdot = zeros(2,1); 
     xdot(1) = - B(2) * x(1) * x(2);    %% d[A]/dt = k[A][C]
     xdot(2) =  - B(1) * x(2) ;   %% d[C]/dt = k'[C]
     dS = xdot;
     end
    S = Sv;
 end
 
%%Least square fitting
B0 = rand(2,1);
[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t(:),[Aconc(:) Cconc(:)], [0;0],[99;99]);
 
%%output parameters value
fprintf(1, '\n')
fprintf(1,'Coefficients:\n')
    for k1 = 1:length(B)
    fprintf(1, '\t\tB(%d) = %e\n', k1, B(k1))
    end
tx = linspace(min(t), max(t));
Plfit = plmodel(B, tx);
 
%%plot graph
figure(1)
plot(t, Aconc/8, 'p');    %plot graph in ratio
hold on
plot(t,Cconc/0.993, 'p'); %plot graph in ratio
hold on
hlp = plot(tx, Plfit(:,1)/8, '-b');
hold on
plot(tx, Plfit(:,2)/0.0993, '-r');
ylim ([0,1])
hold off
xlabel('Time (min)')
ylabel('C/C0')
legend('A','C','fitA','fitC');
end

[Torsten]

Since the maximum values for A are 8 and for C 0.993, lsqcurvefit tries to approximate A "better" than C in order to minimize the sum of squares.

I think this code is more fair in approximating C:

function curveFitting
%Data, t = time; Aconc = concentration of A; Bconc = concentration of B
t = [0 5 10 20 30 45 60];
Aconc = 8*[1;0.546547;0.339339;0.215465;0.175676;0.15991;0.153153]/8;
Cconc = 0.993*[1;0.683173;0.556456;0.448287;0.382121;0.316049;0.257132]/0.993;

    function S = plmodel(B,t) 
        x0 = [8;0.993]; %%Initial concentration
        [T,Sv] = ode45(@DifEq,t,x0); %%Runge Kutta 4th order
        function dS = DifEq(t,x)  
            xdot = zeros(2,1); 
            xdot(1) = - B(2) * x(1) * x(2);    %% d[A]/dt = k[A][C]
            xdot(2) =  - B(1) * x(2) ;   %% d[C]/dt = k'[C]
            dS = xdot;
        end
        S = [Sv(:,1)/8,Sv(:,2)/0.993];
    end

%%Least square fitting
B0 = rand(2,1);
[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t(:),[Aconc(:) Cconc(:)], [0;0],[99;99]);

%%output parameters value
fprintf(1, '\n')
fprintf(1,'Coefficients:\n')
for k1 = 1:length(B)
    fprintf(1, '\t\tB(%d) = %e\n', k1, B(k1))
end
tx = linspace(min(t), max(t));
Plfit = plmodel(B, tx);

%%plot graph
figure(1)
plot(t, Aconc, 'p');    %plot graph in ratio
hold on
plot(t,Cconc, 'p'); %plot graph in ratio
hold on
hlp = plot(tx, Plfit(:,1), '-b');
hold on
plot(tx, Plfit(:,2), '-r');
ylim ([0,1])
hold off
xlabel('Time (min)')
ylabel('C/C0')
legend('A','C','fitA','fitC');
end

[Solomon Thong]

Hi. Thank you for the code. The code works wonders and the explanation really explain the problem that I faced. The fitting seems better if I use Aconc and Bconc in ratio (Aconc/8 and Cconc/0.993) instead of using Aconc and Bconc if use:

        S = [Sv(:,1)/8,Sv(:,2)/0.993];

Thank you very much!

Answer 1

[Star Strider]

If the differential equations are coded correctly, the curves should fit. Check that to be certain that both do, and that they accurately model the process that created the data.

See Parameter Estimation for a System of Differential Equations for a relevant example.

[Solomon Thong]

Hi. Thank you for the answer. It seems like my model does not describe the system accurately. I have one more question regarding the resnorm of the fitting. The resnorm calculate by the system is pretty low (resnorm = 0.0266) despite having one of the curve does not fit the data well. Does it mean that only the best fitted data is used to calculate the resnorm instead of calculating the overall residual of the two different data? Is there a way to know the overall error of the model when more equations are involved when describing the system?

[Star Strider]

My pleasure!

It is possible to get a low ‘resnorm’ result if the magnitudes of the data are also low. The overall error would need to be calculated as the sum of the individual errors. From the lsqcurvefit documentation, it appears that all the errors are combined in that calculation. Using the ‘residual’ output would likely provide the desired informaiton, however I have not specifically examined that.

[Solomon Thong]

Hi. Thank you for the reply. I have changed the order of magnitude of my data (so the order of magnitude is similar) and it seems that the fitting has improve. Thank very much for the explanation!

[Star Strider]

As always, my pleasure!