Trying to fit a differential equation to some data using lsqcurvefit but getting bad results

Hello,

Actually, I am trying to fit a differential equation () to some data using lsqcurvefit. I enetered my initial value in a way that completely fit to the data and I know that it should easily fit because I tried it with cftool. But my code cannot fit the result of the differential equation to my data and gives the following message in command window:

Optimization completed because the size of the gradient at the initial point is less than the value of the optimality tolerance.

Could you please help me to make my differential equation fit to the experimental data?

clc
clear all;
close all;
[d,s,r] = xlsread('Temp.csv');
x = d(:,1);
for i=1:size(x,1)
f1(i)= 0.3148*exp(-((x(i)-67.02)/2.439).^2);
end
f1=f1';
for i=1:size(f1,1)
    %converting temperature to time
    Time(i)=(x(i)-21)/1.5;
end
Time=Time';
B0 = [0.521e6,294.15,1.5];

% [B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@firstorderDSC1,B0,Time,x);
% [B,resnorm,residual] = lsqcurvefit(@firstorderDSC1,B0,Time,f1);
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@firstorderDSC1,B0,Time,f1);
figure;
plot(Time,f1,Time,firstorderDSC1(B,Time))
legend('gaussian','miles')

The following is "firstorderDSC1" function;

function dndt2 = firstorderDSC1(B, t)

[time,n] = ode45(@firstorderDSC, t, 0);
for i=1:size(n,1)
dndt2(i,1) = firstorderDSC(n(i), time(i));
end

    function denaturate = firstorderDSC(n, t)
%           T = (21 + 1.5*t)+273.15;
          a= exp(183.8); %(1/s)
%           delE= 0.521*10^6; %J/mol
%           R=8.3144598; %J/(K*mol)

%           k = a*exp(-delE/(R*T));%+beta*heaviside(T-(55+273.15));

          dndt=a*exp(-B(1)/(B(2)+B(3)*t))*(1-n);

          denaturate=dndt;
      end  
end

I have also attached "Temp.csv"

[Are Mjaavatten]

There are some problems with your problem formulation:

1) Since B occurs only in the combination B(1)/(B(2)+B(3)*t, you cannot determine B(1) and B(2) separately, but only the ratio B(1)/B(2). Therefore you should reduce B to two components: C(1) = B(1)/B(2) and C(2) = B(3).

2) With B = B0, the highest value for the exponent in the expression:

dndt=a*exp(-B(1)/(B(2)+B(3)*t))*(1-n);

is -1437.7 (for t=max(Time)=46) The value of the exponential is then around 1e-769. The smallest absolute value of a double in Matlab is 2.2e-308, thus the exp function returns 0. Your function firstorderDSC1 will then return 0 for a wide range of parameters around B0, so the gradient is 0, and lsqcurvefit will not know in which direction to search.

3) Your f1 array as a function of Time has a peak around Time =30, but is very small elsewhere. Finding parameters that lets firstorderDSC1 mimic this behaviour may not be possible. Are you sure your model is correct?

I recommend that you take a fresh look at your underlying problem and see if you can find an alternative approach.

[Faezeh Manesh]

Hello Are,

Thanks for your comments. I can fix the first problem that you mentioned. but about the last comment, f1 function gradually increases from 0 to 0.3 and then again decreases to 0. In fact f1 is a normal distribution function. I have previously fitted these two together using the following code. But now, I don't know why I cannot use lsqcurvefit to run the same curve-fitting problem:

In fact the following code solves my ODE using ode45 and generates a curve

clc
clear all
tspan = [0 46];
options = odeset('AbsTol', 1e-8,'RelTol',1e-8);%,'OutputFcn',@odeplot);
odefun = @(t,x) firstorderDSC(x,t);

[time,n] = ode45(odefun, tspan, 0, options);
for i=1:size(time,1)
    T(i) = (21 + 1.5.*time(i)); %temperature in Celsius
    dndt2(i)=firstorderDSC(n(i),time(i)); %normalized rate of change of number of molecules
end
T=T';
dndt2=dndt2';
figure;

plot(T,dndt2);
xlabel('Temperature (deg C)');
ylabel('Proportional denaturation rate, dn/dt');

endogram=[time,n];

function denaturate = firstorderDSC(n, t);
%n=n
T = (21 + 1.5*t)+273.15; %T must be in Kelvin here
% n is the proportion denatured (state variable), goes from 0 to 1

beta=1;

a= exp(183.8); %(1/s)
delE= 0.521*10^6; %J/mol
R=8.3144598; %J/(K*mol)

k = a*exp(-delE/(R*T));%+beta*heaviside(T-(55+273.15));


dndt=k*(1-n);
denaturate=dndt;

end

and here I have fitted the result of my ODE with a normal distribution function (1 term gaussian) using cftool and f1 function is this gaussian function that I fitted to the result of my ODE.

In lsqcurvefit, I am trying to do the same thing just with this difference that I have written my ODE in a parametric form (using B(1), B(2) and B(3)) and I set the initial values for these 3 parameters equal to the values that gives me the above results. I don't understand why it doesn't give me the same result.

[Are Mjaavatten]

You write:

I have written my ODE in a parametric form (using B(1), B(2) and B(3)) and I set the initial values for these 3 parameters equal to the values that give me the above results.

Not really:

>> delE= 0.521*10^6;R=8.3144598;T = 294.15;delE/(R*T)
ans =
  213.0271
>> B = [0.521e6,294.15,1.5]; t = 0; B(1)/(B(2)+B(3)*t)
ans =
   1.7712e+03

However, I have a more fundamental question about your problem: Why is temperature T a linear function of time? That does not seem very realistic. My approach would be to model mass and energy balances simultaneously, which would lead to a differential equation of the form:

, where

In Matlab:

z0 = [0;T0];
parameters = something;
fun = @(t,z) reactor(t,z,parameters);
[t,z] = ode45(fun,z0,tspan);

function dxdt = reactor(t,z,parameters)
...
end

Tuneable parameters may for instance be enthalpy of reaction, heat capacities or heating/cooling. You may use lsqcurvefit as before to try to match some measured or expected behaviour.

Good luck!

[Faezeh Manesh]

Hello,

Thanks for your answer and sorry for my late response.

Actually, since I am trying to model DSC test (which is a common test in chemical engineering), in this test temperature is a liner function of time, in fact you choose that how your temperture should change with respect to time.

Anyway, Based on all your feedbacks here I tried some options:

1) I changed my function in a way that only temperature and n are available in moy ODE (so I omitted time from my equation)

2) I decreased my OptimalityTolerance and FunctionTolerance

3) I also tried 0.95*B0 instead of B0 itself

But still I get the same result as before which is a flat line. I don't know what to do to fix the problem. I really appreciate your help.

I have also attached temp.csv

and my code becomes the following

clc
clear all;
close all;
[d,s,r] = xlsread('Temp.csv');

x = d(:,1);
for i=1:size(x,1)
f1(i)=  0.2099*exp(-((x(i)-340.2)/2.439).^2);
end
f1=f1';

B0 =[-5.9529e+04];

options = optimoptions('lsqcurvefit','OptimalityTolerance',1e-12,'FunctionTolerance',1e-9);
lb = [];
ub = [];
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@firstorderDSC1,B0,x,f1,lb,ub,options);

figure;
plot(x,f1,x,firstorderDSC1(B,x))
legend('gaussian','miles')

and my firstorderDSC1 function as:

function dndt2 = firstorderDSC1(B, T)

options = odeset('AbsTol', 1e-8,'RelTol',1e-8);
[temp,n] = ode45(@firstorderDSC, T, 0, options);
for i=1:size(n,1)
dndt2(i,1) = firstorderDSC(n(i), temp(i));
end

    function denaturate = firstorderDSC(n, T)

          a= (exp(183.8)/1.5); 

          dndT=a*exp(B(1)/T)*(1-n);
          denaturate=dndT;
      end  
end

Answer 1

[Alan Weiss]

My interpretation of your reported lsqcurvefit result is that the solver agrees that you started with the optimal point. It didn't take a step because there was nowhere better to go.

To check whether I am correct, start lsqcurvefit from a different point than the optimal one and see if it takes you back to the optimal point. For example, try

B0 = 0.95*B0;

before you call lsqcurvefit.

Alan Weiss

MATLAB mathematical toolbox documentation

[Faezeh Manesh]

Hello Alan,

Thanks for your response and sorry for my late reply.

I did that and the result didn't change.

I explained it in detail in last comment.

If you have any other idea, i really appreciate it.

regards,

Faezeh