Parameter fit differential equation to dataset
33 次查看(过去 30 天)
显示 更早的评论
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);
0 个评论
采纳的回答
Star Strider
2019-10-10
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 个评论
Star Strider
2019-10-11
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
2019-10-11
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.
更多回答(1 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Interpolation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!