Parameter fit differential equation to dataset

30 vues (au cours des 30 derniers jours)
Maria Alvarez
Maria Alvarez le 10 Oct 2019
Commenté : Star Strider le 11 Oct 2019
Hello!
I am trying to fit a differential equation to a dataset by adjusting the equation's parameters (c(1) and c(2)). My datasets consist of an evolution of sealevel (s) over time (t) and temperature (temp) over time. The equation is as follows:
ds/ dt=-(1-c(1))*c(2).*(tempi(t))
My code is not working and was hoping someone could please give me some pointers.
CODE:
function bestc = odeparamAsmb5()
%Initial data
load rcp85_antsmbmidmatlab.txt;
load rcp85_temperaturemidmatlab.txt;
time= rcp85_temperaturemidmatlab(:,1)-2006+1;%modifying the vector from [2006 2100] to [1 95] to match "i"
temp= rcp85_temperaturemidmatlab(:,2)+0.8;
sealevel=rcp85_antsmbmidmatlab(:,2);
plot(time,sealevel)
%information
tempi=@(t) interp1(time,temp,t);
f=@(c,t)-(1-c(1))*c(2).*(tempi(t));
simsealevel= zeros(1,95);
for t=1:1:95
simsealevel(t)= integral(f,t,t+1);
hold on
plot(time,simsealevel, '-r')
end
% Optimise difference between simsealevel and sealevel
fun= @(c,t)-(1-c(1))*c(2).*(tempi(t))- sealevel;
c0= [0.25,0.15];
lb= [0.25;0.15];
ub= [0.25;0.35];
x = lsqnonlin(fun,c0,lb,ub);
ERRORS:
Not enough input arguments.
Error in odeparamAsmb5>@(c,t)-(1-c(1))*c(2).*(tempi(t)) (line 14)
f=@(c,t)-(1-c(1))*c(2).*(tempi(t));
Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);
Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);
Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);
Error in odeparamAsmb5 (line 20)
simsealevel(t)= integral(f,t,t+1);

Réponse acceptée

Star Strider
Star Strider le 10 Oct 2019
I have discussed this several times over the years. See for example Monod kinetics and curve fitting and Parameter Estimation for a System Of Differential Equations. I will be glad to help you, however I will need your data.
I suggest that you do not interpolate your data. The reason is that it creates data points that may not actually exist, since you have no idea what your data are doing except at the times and values you measured them.
  6 commentaires
Star Strider
Star Strider le 11 Oct 2019
Maria Alvarez’s Answer moved here —
Thank you, Star Strider!
Yes, c(1) and c(2) are the parameters that I am trying to estimate (each with its lower and upper bounds). However, "tempi (t)" is "T(t)", Temperature as a function of time. The issue is that I don't have a mathematical representation for T(t), just the dataset of Temperature over time.
Star Strider
Star Strider le 11 Oct 2019
This runs, however it will need some help from you to make the fit better. What you are doing is far outside my areas of expertise, so I cannot help you with the intracicies of the objective function.
Try this:
rcp85_antsmbmidmatlab = load('rcp85_antsmbmidmatlab.txt');
rcp85_temperaturemidmatlab = load('rcp85_temperaturemidmatlab.txt');
time= rcp85_temperaturemidmatlab(:,1)-2006+1;%modifying the vector from [2006 2100] to [1 95] to match "i"
temp= rcp85_temperaturemidmatlab(:,2)+0.8;
sealevel=rcp85_antsmbmidmatlab(:,2);
function S = objfun(c,t,temp)
s0 = c(3);
c = [c(1), c(2)];
[T,S]=ode45(@DifEq,t,s0);
function dsdt = DifEq(t,s)
tempi = @(t) interp1(time,temp,t);
dsdt = -(1-c(1))*c(2).*tempi(t);
end
end
c0= [0.255,0.155,0.1];
lb= [0.25;0.15];
ub= [0.26;0.35];
[cest,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@(c,t)objfun(c,t,temp),c0,time,sealevel,lb,ub);
seaest = objfun(cest,time,temp);
figure
plot3(time, temp, sealevel)
hold on
plot3(time, temp, seaest)
hold off
grid on
xlabel('Time')
ylabel('Temp')
zlabel('Sea Level')
The problem with the fit is likely due to the parameter constraints, since I can get a very good fit with no constraints at all.
The unconstrained parameter estimates are:
cest =
0.996983010442519 0.056715050093189 0.002671958952739
for ‘c(1)’, ‘c(2)’ and the initial condition ‘s0’, respectively, with the norm of the residuals:
Rsdnrm =
1.083258730496945e-04
So the fit is quite good.
The constrained fit is much worse, and the estimate actually goes in the opposite direction from the observed data with respect to both temperature and time.
Experiment to get the result you want.

Connectez-vous pour commenter.

Plus de réponses (1)

Maria Alvarez
Maria Alvarez le 11 Oct 2019
Thanks again Star Strider!!

Community Treasure Hunt

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

Start Hunting!

Translated by