How to fit data with a solution of a non linear differential equation containing multiple unknown coefficient

Hi everybody,

So it's been quite some time since I have used Matlab and now I have had a hard time figurint out how to proceed with my problem.

So I have measured some data and I would like to fit it with a certain function.

This latter function is found by solving numerically a non linear differential equation with 6 unknown coefficients.

What I want to do is by using a least-squares method (lsqcurvefit) I would vary the 6 unknown coefficients of the non linear ODE so that after the loop routine the solution to this equation fits my data.

And I'm stuck in thinking how to label all this. I understand the examples of the lsqcurvefit but it applies to a function. And I want to apply it to a equation so that the solution fits the data.

Here is the equation with the unknown coefficients a, b1, b2, c1, c2 and d. I aleady have a guess for all these values and a range in which they should lie within.

When I solve this equation numerically in matlab (with certain values), I can plot the function u(t). But what I want to plot and fit is actually .

So in summary: I want to write a routine loop varying my 6 coefficients in my non linear equation in order for the p(t) function to fit my data.

Any suggestion from you guys ?

Thanks a lot.

Best,

Thibaut

Answer 1

[Torsten]

In cases like yours, we always refer to

https://de.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

Best wishes

Torsten.

[Thibarely]

Thanks for the help.

I started with what you gave but so far was not successful as my Matlab skills are still limited :)

Here is what I've written:

%calling statement
Time=[0,0.0625 0.125 0.1875 0.25 0.3125 0.3750 0.4375 0.5 0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1 1.0625 1.125 1.1875 1.25 1.3125 1.375 1.4375 1.5 1.625 1.75 1.875]; 
PLdata=[1 0.8 0.6 0.57 0.53 0.5 0.45 0.4 0.35 0.3 0.27 0.23 0.2 0.18 0.16 0.13 0.115 0.08 0.07 0.07 0.06 0.05 0.05 0.06 0.06 0.045 0.04 0.035];
B0=rand(6,1)*100; %creating a vector for the 6 unknown parameters
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@plmodel,B0,Time,PLdata); %initiating the lsqcurvefit loop

Then I created my plmodel.m in the same folder:

function U = plmodel(B,t)
%variables: x=U pl=P
%Parameters: p_0=B(1), tau_p=B(2), tau_n=B(3), c_n=B(4), c_p=B(5), k*=B(6)

x0=rand(2,1);
[t,UP]=ode45(@DifEq,t,x0);

    function dU=DifEq(t,x)
        xdot=zeros(1,1);
        pl=zeros(1,1);
        xdot=-x.*(x+B(1)).*(1/(B(2).*x+B(3).*(x+B(1)))-B(4).*x-B(5).*(B(1)+x)-B(6));
        pl=B(6).*x.*(B(1)+x);
        dU=xdot;
    end

U=UP(:,1);
end

I have to say, I am a bit lost here. Anything more specific to help me understand how this loop works ?

Thanks !

Answer 2

[Torsten]

function main
  Time = [0,0.0625 0.125 0.1875 0.25 0.3125 0.3750 0.4375 0.5 0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1 1.0625 1.125 1.1875 1.25 1.3125 1.375 1.4375 1.5 1.625 1.75 1.875]; 
  PLdata = [1 0.8 0.6 0.57 0.53 0.5 0.45 0.4 0.35 0.3 0.27 0.23 0.2 0.18 0.16 0.13 0.115 0.08 0.07 0.07 0.06 0.05 0.05 0.06 0.06 0.045 0.04 0.035];
  B0 = rand(6,1)*100; %creating a vector for the 6 unknown parameters
  [B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@plmodel,B0,Time,PLdata); %initiating the lsqcurvefit loop
end

function u = plmodel(B,t)
  %variables: x=U pl=P
  %Parameters: p_0=B(1), tau_p=B(2), tau_n=B(3), c_n=B(4), c_p=B(5), k*=B(6)
  u0 = 1.0;   %you must specify the correct starting value u0 for u !!!
  [~,u] = ode45(@(tf,uf)DifEq(tf,uf,B),t,u0);
  u = B(6)*u.*(B(1)+u);
end

function du = DifEq(t,u,B)
  du = -u*(-u+B(1))*(1/(B(2)*u+B(3)*(u+B(1)))-B(4)*u-B(5)*(B(1)+u)-B(6));
end

[Thibarely]

Hello Torsten. Thank you again for the help. I understand how I should include my "u" function in the whole function but I do not understand all of your code so I decided to follow again these topics: https://fr.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting and https://fr.mathworks.com/matlabcentral/answers/314965-parameter-estimation-for-a-system-of-differential-equations#comment_411192

Here is what I've came up with so far:

function trplfit

%This function will fit the trpl curves in order to obtain the 6 unknown
%coefficients 

function u= plmodel(B,t)
   %Parameters: p_0=B(1), tau_p=B(2), tau_n=B(3), c_n=B(4), c_p=B(5), k*=B(6)
  u0 = [1];
  [T,Uv] = ode45(@DifEq,t,u0);
  u= B(6).*u.*(B(1)+u);
  
  function dU = DifEq(t,u)
  udot = zeros(1,1);    
  udot(1) = -u(1).*(-u(1)+B(1)).*(1/(B(2).*u(1)+B(3).*(u(1)+B(1)))-B(4).*u(1)-B(5).*(B(1)+u(1))-B(6));
  dU=udot;
  end
U=Uv;    
end

t = [0 0.0625 0.125 0.1875 0.25 0.3125 0.3750 0.4375 0.5 0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1 1.0625 1.125 1.1875 1.25 1.3125 1.375 1.4375 1.5 1.625 1.75 1.875]; 
u = [1 0.8 0.6 0.57 0.53 0.5 0.45 0.4 0.35 0.3 0.27 0.23 0.2 0.18 0.16 0.13 0.115 0.08 0.07 0.07 0.06 0.05 0.05 0.06 0.06 0.045 0.04 0.035 0.032];
B0 = [3*10^15 871*10^-9 511*10^-9 4*10^-29 4*10^-29 4.78*10^-11]; %creating a vector for the 6 unknown parameters
[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t,u); %initiating the lsqcurvefit loop
end

But I get these errors:

Undefined function or variable 'u'.

Error in trplfit/plmodel (line 10)
  u= B(6).*u.*(B(1)+u);

Error in lsqcurvefit (line 202)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in trplfit (line 23)
[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t,u); %initiating the lsqcurvefit loop

Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

What would you change in the code ?

Thanks

[Torsten]

[T,u] = ode45(@DifEq,t,u0);

and delete line

U=Uv

[Thibarely]

Thank you again.

I've proceeded with your suggestion. (also I removed one number in the u vector to match the t vector size)

I now have the following errors:

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

Error in trplfit (line 26)
[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t,u); %initiating the lsqcurvefit loop

How is my function represented ? Because I have the feeling that I have on one side a vector (ydata) and on the other side I have something different and therefore the calculation for lsqcurvefit cannot start.

What would you recommend ?

Best,

Thibaut

[Torsten]

In the call

[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t,u); 

t and u must have the same length.

In plmodel, t and u also must have the same length.

And note that you start with u0=1. So t=0 should be included in your t-vector.

Answer 3

[Thibarely]

The script works now but does not take into account that I want to use the least square method to fit

which in my code is given by:

u= B(4).*u.*(3*10^15+u);

So what the script does is that it fits the solution of the non linear equation and then uses the values found inside

u= B(4).*u.*(3*10^15+u);

Therefore the plot I get is completely wrong.

How can I change my code so that the lsqcurvefit is applied on the last function and not on the solution of the equation ?

Here is my code:

function trplfit

%This function will fit the trpl curves in order to obtain the 6 unknown
%coefficients 

function u=plmodel(B,t)
   %Parameters: p_0=B(1), tau_p=B(2), tau_n=B(3), c_n=B(4), c_p=B(5), k*=B(6)
u0 = 1;
[T,u] = ode45(@DifEq,t,u0);

  
  function dU = DifEq(t,u)
  udot = 1;    
  udot = -u.*(u+3*10^15).*(1/(B(1).*u+B(2).*(u+3*10^15))-B(3).*u-B(4));
  dU=udot;
    end
  u= B(4).*u.*(3*10^15+u);  
end


t = 1*10^-6*[0 0.0625 0.125 0.1875 0.25 0.3125 0.3750 0.4375 0.5 0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1 1.0625 1.125 1.1875 1.25 1.3125 1.375 1.4375 1.5 1.625 1.75 1.875] 
u = [1; 0.8; 0.7; 0.6; 0.55; 0.5; 0.45; 0.4; 0.35; 0.3; 0.27; 0.23; 0.2; 0.18; 0.16; 0.13; 0.115; 0.118; 0.09; 0.092; 0.068; 0.07; 0.045; 0.045; 0.047; 0.03; 0.026; 0.022]
B0 = [871*10^-9; 511*10^-9; 8.83*10^-29; 4.78*10^-11] %creating a vector for the 4 unknown parameters
plmodel(B0,t)

[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t,u); %initiating the lsqcurvefit loop



fprintf(1,'\tCoefficients:\n')
for k1 = 1:length(B)
    fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
end

tv = linspace(min(t), max(t)); 
PLfit = plmodel(B, tv);

figure(1)
semilogy(t, u, 'p');
hold on
hlp = semilogy(tv, PLfit);
hold off
xlabel('Decay time (ns)')
ylabel('Normalized Photoluminescence')
legend(hlp, 'fit')

end

[Torsten]

I don't get what you mean.

What you do in your code is to fit your input vector "u" to the output vector "u" of "plModel", thus to B(4).*udiff.*(3*10^15+udiff) where "udiff" is the solution of the differential equation

d udiff/dt = -udiff.*(-udiff+3*10^15).*(1/(B(1).*udiff+B(2).*(udiff+3*10^15))-B(3).*udiff-B(4)) 

If you want to fit your input vector u directly to udiff obtained from the solution of the differential equation, remove the line

u= B(4).*u.*(3*10^15+u);  

in your code.

And please don't create new "Answers" every time you make "Comments".

[Thibarely]

My bad I did it only once but you do seem to be a little grumpy. That's ok.

Here is what I meant:

I can see that my lsqcurvefit is fitting my data using the solution of the non linear equation.

udot = -u.*(u+3*10^15).*(1/(B(1).*u+B(2).*(u+3*10^15))-B(3).*u-B(4));

However as I said, I wanted the lsqcurvefit to fit my data using a formula expressed using the solution of this non linear equation:

Say "u" is the solution of my non linear equation. I want the lsqcurvefit to fit my data not with "u" but with

B(4).*u.*(3*10^15+u)

Is it clearer now ?

Thanks again

[Torsten]

Say "u" is the solution of my non linear equation. I want the lsqcurvefit to fit my data not with "u" but with

B(4).*u.*(3*10^15+u)

But that's exactly what lsqcurvefit does in your code:

If u is the solution of the differential equation

udot = -u.*(-u+3*10^15).*(1/(B(1).*u+B(2).*(u+3*10^15))-B(3).*u-B(4));

, lsqcurvefit fits B(4)*u*(3e15+u) to your data vector

[1; 0.8; 0.7; 0.6; 0.55; 0.5; 0.45; 0.4; 0.35; 0.3; 0.27; 0.23; 0.2; 0.18; 0.16; 0.13; 0.115; 0.118; 0.09; 0.092; 0.068; 0.07; 0.045; 0.045; 0.047; 0.03; 0.026; 0.022]

[Torsten]

Although then I don't understand t(1)=0, u(1)=1 in your data vector. If you wanted to fit against B(4)*u*(3e15+u), your first u(1) should be much larger than 1. Seems we get our wires crossed.

[Thibarely]

Hi Torsten,

Yes I also thought of that, and I changed my u0 to something more reasonable.

Also, the formula that I wrote down was wrong (two negative signs instead of two positive signs) so now I can get a fit but this fit does not match my data and I think I know why but I need your help to see how I can make it work.

Here is my code with my data:

function trplfitmapi39
options=optimoptions(@lsqcurvefit, 'MaxFunctionEvaluations',2000000, 'MaxIterations', 200000)
format shortg
%This function will fit the trpl curves in order to obtain the 6 unknown
%coefficients 

function u=plmodel(B,t)
   %Parameters: tau_p=B(1), tau_n=B(2), C=B(3), k*=B(4)
u0 = 6*10^13;
[T,u] = ode45(@DifEq,t,u0);
  
  function dU = DifEq(t,u)
  udot = 0;    
  udot = -u.*(u+B(1)).*(1/(B(2).*u+B(3).*(u+B(1)))+B(4).*u+B(5));
  dU=udot;
  end
  
u= B(5).*u.*(B(1)+u);

end

load mapi39b.csv;
t = mapi39b(:,1);
u = mapi39b(:,2);
B0 = [1*10^13; 871*10^-9; 511*10^-9; 8*10^-29; 4*10^-11] %creating a vector for the unknown parameters

c=plmodel(B0,t)
figure(1)
semilogy(t,c);
%lb= [1*10^12; 1*10^-9; 1*10^-9; 8*10^-31; 4*10^-13]
%ub= [1*10^16; 1*10^-5; 1*10^-5; 8*10^-26; 4*10^-9] 

[B,resnorm,residual,exitflag,output] = lsqcurvefit(@plmodel,B0,t,u,1,+inf,options); %initiating the lsqcurvefit loop


fprintf(1,'\tCoefficients:\n')
for k1 = 1:length(B)
    fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
end

tv = linspace(min(t), max(t)); 
PLfit = plmodel(B, t);

figure(2)
semilogy(t, u, 'o');
hold on
hlp = semilogy(t, PLfit);
hold off
xlabel('Decay time (ns)')
ylabel('Normalized Photoluminescence')
legend(hlp, 'fit')

end

I think that in order to get a correct fit I have to normalize both my data and the function "u" so that the lsqcurvefit will be able to work.

How do I do that ? I mean normalizing a vector for my data seems ok, and I can even do it before. However I've tried normalizing the function "u" with normalize(u) but this does not work as it is a function.