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

13 views (last 30 days)
Faezeh Manesh on 24 Feb 2020
Commented: Faezeh Manesh on 4 Mar 2020
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.
clc
clear all;
close all;
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"

Show 1 older comment
Faezeh Manesh on 25 Feb 2020
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 on 27 Feb 2020
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 on 4 Mar 2020
Hello,
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;
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

Alan Weiss on 25 Feb 2020
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

#### 1 Comment

Faezeh Manesh on 4 Mar 2020
Hello Alan,