Parameter fit differential equation to dataset

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);

Answer 1

[Star Strider]

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.

[Maria Alvarez]

Thanks very much, Star Strider!! I will check those examples out and try to sort them out myself. Duly noted on the data interpolation

[Star Strider]

As always, my pleasure!

[Maria Alvarez]

Hi again Star Strider,

I went through the Mond Kinetics example and am not sure whether using ODE45 would apply to my current differential equation (as the ds/dt= -(1- c1) * c2 * T(t) equation does not include "s" as a variable) and was hoping you could please clarify.

I had modified a code to fit a different equation (ds/dt= (c1* T(t) - s(t) ) / c2) and dataset and it worked (please see CODE 1 below and its datasets). Initially, I tried modifying this code to fit the new differential equation (which no longer includes "s" as a variable) to a different dataset but the fit was very poor (please see CODE 2 below, also attaching the datasets) which is why I decided to go a different way around it with CODE 3 (please see below).

Any pointers on how to go around this would be greatly appreciated!!

CODE 1:

function bestc = odeparamThexp()
%Initial data

load rcp85_expansionmidmatlab.txt;
load rcp85_temperaturemidmatlab.txt;

time= rcp85_temperaturemidmatlab(:,1)-2006+1;%modifying the vector from [2006 2100] to [1 95] to match tSpan
temp= rcp85_temperaturemidmatlab(:,2)+0.8;
sealevel=rcp85_expansionmidmatlab(:,2);
plot(time,sealevel)

%ODE information
tSpan = [1 95];
z0 = [0.0118];
tempi=@(t) interp1(time,temp,t);

%Initial guess

c0= [0.2,82];

fun=@(t,z)updateStates(t,z,c0,tempi);
ODE_Sol= ode45(fun,tSpan,z0); % Run the ODE
simsealevel = deval(ODE_Sol, time); % Evaluate the solution at the experimental time steps

hold on
plot(time, simsealevel, '-r')

%% Set up optimization
myObjective = @(x) objFcn(x,time,sealevel,tempi,tSpan,z0);
lb = [0.2,81.7];
ub = [0.626,1290];

bestc = lsqnonlin(myObjective,c0, lb, ub);


%% Plot best result
ODE_Sol = ode45(@(t,z)updateStates(t,z,bestc,tempi),tSpan, z0);

bestsealevel = deval(ODE_Sol, time);

plot(time, bestsealevel, '-g')
title('Thermal expansion over time')
xlabel('t') 
ylabel('SLR [m]') 
legend('IPCC Data','Initial Param','Best Param');


function f = updateStates(t,z,c,tempi)


f = (c(1).*tempi(t)-z)*(1/(c(2)));


function cost= objFcn (x,time,sealevel,tempi,tSpan,z0)


ODE_Sol = ode45(@(t,z)updateStates(t,z,x,tempi), tSpan, z0);
simsealevel = deval(ODE_Sol, time);


cost = simsealevel-sealevel;

CODE 2:

function bestc = odeparamAsmb()
%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 tSpan
temp= rcp85_temperaturemidmatlab(:,2)+0.8;
sealevel=rcp85_antsmbmidmatlab(:,2);
plot(time,sealevel)

%information
tSpan = [1 95];
z0 = [-0.0001];
tempi=@(t) interp1(time,temp,t);

%Initial guess

c0= [0.25,0.0];

fun=@(t,z)updateStates(t,z,c0,tempi);
ODE_Sol= ode45(fun,tSpan,z0); % Run the ODE
simsealevel = deval(ODE_Sol, time); % Evaluate the solution at the experimental time steps

hold on
plot(time, simsealevel, '-r')

%% Set up optimization
myObjective = @(x) objFcn(x,time,sealevel,tempi,tSpan,z0);
lb = [0.25,0.15];
ub = [0.25,0.35];

bestc = lsqnonlin(myObjective,c0, lb, ub);


%% Plot best result
ODE_Sol = ode45(@(t,z)updateStates(t,z,bestc,tempi),tSpan, z0);

bestsealevel = deval(ODE_Sol, time);

plot(time, bestsealevel, '-g')
legend('IPCC Data','Initial Param','Best Param');


function f = updateStates(t,z,c,tempi)


f = -(1-(c(1))*c(2).*(tempi(t)));


function cost= objFcn (x,time,sealevel,tempi,tSpan,z0)


ODE_Sol = ode45(@(t,z)updateStates(t,z,x,tempi), tSpan, z0);
simsealevel = deval(ODE_Sol, time);


cost = simsealevel-sealevel;

CODE 3:

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);

[Star Strider]

Your original differential equation:

ds/dt = -(1-c(1))*c(2).*(tempi(t))

simplified to:

ds/dt = -(1-c(1))*c(2).*t

would integrate to :

s(t) = (c(2)*(c(1) - 1)*t^2)/2 + s(0)

with ‘s(0)’ being the initial condition, a parameter you can easily include in the estimation.

This appears to be straightforward, and no other integrations are necessary. Are ‘c(1)’ and ‘c(2)’ the parameters you are estimating (in addition to ‘s(0)’), with your independent variable being ‘t’?

[Star Strider]

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]

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.

Answer 2

[Maria Alvarez]

Thanks again Star Strider!!

[Star Strider]

As always, my pleasure!