ode45 + lsqcurvefit- multiple ODEs

Hi,

I am trying to fit my experimental data to a model I have created that include 2 coupled ODEs. I created a separate function for defining the 2 ODEs, where a is the array that contains the 2 constants I'm trying to fit:

function dy = CalcCoupODEs(t,y,a) %dy= matrix with all 2 ODEs

dy = zeros(2,1); % variable definition

dy(1) = 1.41.*a(1) .* ((y(2)*y(1)^2)^2)*y(1) - a(2) .* y(1)*t ; %ODE for R

dy(2) = -a(1)*(y(2)*y(1)^2)^2; %ODE for n

end

Then, I call the function in the following way:

cycle_DMACl = [20, 25, 30, 40, 45];  % this is t variable

R_DMACl = [18.34, 11.29, 7.09, 6.51, 4.396]; %this is y(1)

n_DMACl = [92.31, 61.9, 53.82, 26.04, 11.87]; % this is y(2)

a0 = [1, 1];  %initial guess for constants

lb1 = [0, 0]; %lower bound

a = lsqcurvefit(@CalcCoupODEs, a0, cycle_DMACl(:), R_DMACl(:), n_DMACl(:), lb1);

I get the following error: Warning: Length of lower bounds is > length(x); ignoring extra bounds.

It's obviously not trying to fit the correct thing (a and lb are both 1x2)

What am I doing wrong and what is the b est way to do this?

Thank you!

Answer 1

[Star Strider]

Similar to Monod kinetics and curve fitting (link), I formulated your problem as:

cycle_DMACl = [20, 25, 30, 40, 45];  % this is t variable
R_DMACl = [18.34, 11.29, 7.09, 6.51, 4.396]; %this is y(1)
n_DMACl = [92.31, 61.9, 53.82, 26.04, 11.87]; % this is y(2)
Rn_DMACl = [R_DMACl(:), n_DMACl(:)];                                % Create Single Matrix

a0 = [1, 1];  %initial guess for constants
lb1 = [0, 0]; %lower bound

a = lsqcurvefit(@(b,t)objfcn(b,t), [a0 0 0], cycle_DMACl(:), Rn_DMACl) 

      function yest = objfcn(b,t)
          a(1) = b(1);
          a(2) = b(2);
          ic = [b(3); b(4)];

          function dy = CalcCoupODEs(t,y) %dy= matrix with all 2 ODEs
              dy = zeros(2,1); % variable definition
              dy(1) = 1.41.*a(1) .* ((y(2).*y(1).^2).^2).*y(1) - a(2) .* y(1).*t; %ODE for R
              dy(2) = -a(1).*(y(2).*y(1).^2).^2; %ODE for n
          end

          [t,yest] = ode15s(@CalcCoupODEs, t, ic);

      end

Rn_fit = objfcn(a,cycle_DMACl(:))

figure(1)
plot(cycle_DMACl(:), Rn_DMACl, 'p')
hold on
plot(cycle_DMACl(:), Rn_fit, '-r')
hold off
grid

The problem is that your ‘CalcCoupODEs’ keeps crashing, with the error:

Warning: Failure at t=2.000000e+01.  Unable to meet integration tolerances without reducing the step size below
the smallest value allowed (5.684342e-14) at time t.

So it is encountering a singularity, and the integration stops. You need to fix that. I do not know what system of differential equations you are coding (you did not post them).

My code should work otherwise, although I cannot be certain of that, since the ODE integration crashes on the first time step, even with non-zero and non-identical initial conditions.