MATLAB Answers

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

13 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

Show 1 older comment
Faezeh Manesh
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
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
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.

Sign in to answer this question.


Translated by