Fitting data with ODE equation and multiple inputs

Hello,

I have an equation that needs to be fitted for multiple parameters but am genuinely stuck on how to do it. Please see the equation below:

It is important to note that z is a function of x and I have data for z(x) and dy/dx.

I am interested in fitting parameters for a, b, and n, but the problem does not seem entirely straightforward to me in MATLAB and I am hoping to get some guidance. I've looked at ways to try to fit the ODE using the guide written here and even tried looking at some analagous problems with Monod kinetics, but the way that the solutions have been posed for Monod kinetics is not entirely clear to me since they have data for all of the variables.

Below is some example data that can be used to try and fit.

x = [0 0.0258 0.0515 0.0773 0.1031 0.1289 0.1546]

dy/dx = [3.4301 5.5175 6.3542 5.9365 4.8554 3.6837 2.7136]

z = [943.0802 1136.4406 1308.0878 1439.5741 1533.7588 1602.7919 1648.6146]

Thank you in advance!

[Star Strider]

‘... z is a function of x ...’

Is it possible for us to know what that functional relationship is?

[Dylan Kline]

Well, the data is experimental so x and z are measured independently, though an equation could be fit if that is what you are referring to?

Answer 1

[Star Strider]

It turns out that your differential equation has an analytical solution:

syms a b n x y(x) z y0 
Eq = diff(y) == a*(1-y)^n * exp(-b/z);
ys = dsolve(Eq, y(0) == y0)
ys = simplify(ys, 'Steps',500)
yfcn = matlabFunction(ys, 'Vars',{[a,n,b,y0],x,z})

In the 'Vars' name-value pair, this [a,n,b,y0] concatenates those variables (parameters here) into a single row vector (in that order) called ‘in1’. That is the parameter vector I later refer to as ‘p’.

Using ‘yfcn’ as calculated:

x = [0	0.0258	0.0515	0.0773	0.1031	0.1289	0.1546];
dydx = [3.4301	5.5175	6.3542	5.9365	4.8554	3.6837	2.7136];
z = [943.0802	1136.4406	1308.0878	1439.5741	1533.7588	1602.7919	1648.6146];

xz = [x(:) z(:)];

yfcn = @(in1,x,z)-((in1(:,2)-1.0).*((-in1(:,4)+1.0).^(-in1(:,2)+1.0)./(in1(:,2)-1.0)+in1(:,1).*x.*exp(-in1(:,3)./z))).^(-1.0./(in1(:,2)-1.0))+1.0;

p0 = rand(1,4);

dydxfcn = @(p,xz) gradient(yfcn(p,xz(:,1),xz(:,2))./gradient(xz(:,1)));
B = lsqcurvefit(dydxfcn, p0, xz, dydx(:));

If you want to integrate ‘dydx’, the code changes to:

yv = cumtrapz(x, dydx);

objfcn = @(p,xz) yfcn(p,xz(:,1),xz(:,2));
B = lsqcurvefit(objfcn, p0, xz, yv(:));

with no other changes.

Both of these run without error. I leave it to you to determine if it/they produce the correct results. Note that nonlinear parameter estimation problems can be very sensitive to the initial parameter estimates, and may end up in a local minimum rather than the desired global minimum, with unacceptable parameters. There are ways to deal with this (the ga function being my favourite because I’ve had the most experience with it over the years), however even those do not always find the best parameter set.

[Dylan Kline]

I must say, this is an elegant solution and I genuinely appreciate your time in helping me with this. Also, thank you for being such an active member of this community... I saw your name on almost all of the posts with similar questions and you have truly been a help to us at large.

[Star Strider]

As always, my pleasure!

I very much appreciate your compliment!