Estimating initial conditions of an ODE system using least squares

Question

I am trying to estimate the parameter values of the differential equation,
${dN\over dt}=\lambda N$. I know how to estimate the parameter value of $\lambda$ but I am having trouble in estimating the initial value of $N$, $N_0$. I am using non-linear least squares to estimate the parameters using Matlab through the function lsqnolin. The code is as below and I would like to know if the way I am estimating the initial condition is correct. The actual model is more complex and the data is different but I want to clarify of a way to estimate ODE initial conditions. With the below code I am not that sure if the method is right as if I start with an initial guess of $10^{12}$ for $N_0$ the estimated output is also the same, so I am not sure if is being optimised properly.

function [estimates,error]= leastSquaresModel()

t=[0;1;2;3];

dataN=[10;9;8;7];

initial=[10,10^12];%initial guesses of the estimates and intial condition
options=optimoptions(@lsqnonlin,'MaxFunctionEvaluations',10000,'MaxIterations',10000);
[estimates,error]=lsqnonlin(@(xEstimate)errorFun(xEstimate,dataN,t),initial,[],[],options);

end

function err=errorFun(xEstimate,y_meas,thours)
    y_est=model2(xEstimate,thours);%y_est are the estimated results from the model fit
    err=y_meas-log10(y_est);%y_meass are the data
end


function output= model2(xEstimate,thours)
    intVals=xEstimate(2);  %inital value of ODE (N_0)  
    options = odeset('AbsTol',1e-13,'RelTol',1e-11);

    [~,values] = ode45(@(t,y)Equations(t,y),thours,intVals,options);

    function s=Equations(~,y)

        r=xEstimate(1);
        s=zeros(1,1);
        s(1)=r*y(1);
    end

output=values';

end

Answer 1

The answer consists of two parts. The first the Matlab code error, and the second is an optimization trick.

1. The line [~,values] = ode45(@(t,y)Equations(t,y),thours,intVals,options); gives you a vector. Then you transpose it as output=values';, and you have a row. Thus, when you compute err=y_meas-log10(y_est);, you have a matrix. It is not probably what you expect since what you need is a vector of errors. How to fix it? First, you can replace output=values'; with output=values;. Second, you can control your dimensions using size, isrow, and iscolumn.

2. Your initial guess for the second parameter is $10^{12}$, that is way larger than the initial guess for the first one, that is $10$. It can be a good idea to redefine as initial=[10,12]; and update the model2 function as intVals=10^xEstimate(2);. I.e., it can be nice to use as the second parameter not the value itself, but the power.

With these two modifications, the optimization result is $\lambda=-2.306$ and $N_0=10^{10}$, that is very close to the real value $\lambda=\log(0.1)=-2.303$.

And one more remark. For your model, it is reasonable to assume that $\lambda<0$, i.e., $N$ does not diverge. Then it is a good idea to take the initial guess negative as well.