MATLAB Answers

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

3 views (last 30 days)
Faezeh Manesh
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.
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"
  4 Comments
Faezeh Manesh
Faezeh Manesh on 4 Mar 2020
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

Sign in to comment.

Answers (1)

Alan Weiss
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
Faezeh Manesh on 4 Mar 2020
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

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by