Finding best parametric function estimation for ODE of first order

Hello to every one.

Until now I've solved differential equations numerically under given formulae for their coefficients and given initial / boundary conditions. But recently I've been dealt with the inverse problem, namely: find the best fit to the solution of the ODE regarding some unknown paramateres. The hard part here is due to the fact that we need to determine not a single value but the values of a whole parametric function.

Consider the following initial value problem for ODE of first order:

Solving the equation by separation of the variables we obtain the exact solution

Let's consider however this problem statement: given the Cauchy's problem

whose solution is known over the grid Find the estimated values of the function fitting the solution , i.e. establish that they're eventually located close enough to the curve .

I read a lot about the topic but I didn't find what I searched for. Could you give me some short explanations or guidelines on how to perform this idea in general and later using MATLAB toolboxes?

Kind regards,

Lyudmil Yovkov, PhD

Answer 1

[Torsten]

Search for a book on optimal control of ODEs and DAEs, e.g.

Matthias Gerdts:

Optimal control of ODEs and DAEs

De Gruyter 2011

Or start with an introductory script:

silo.tips/download/optimal-control-of-odes-introductory-lecture

Finally, you will have to use Matlab's lsqcurvefit or fmincon to determine the coefficients c_i on the grid x_i that minimize

sum_{i=0}^{i=n} (u_i - u_i(c_i))^2

where u_i(c_i) are the values for u obtained from the ODE by using the c_i as coefficients in front of du/dx.

[Lyudmil Yovkov]

@Torsten, I shall definitely try what you're suggesting me. Could you give me any direct links to how I may use the optimization procedures combined with ODE solvers? Actually I solve the equations (ODEs and PDEs) using Finite Difference Method and implementing the famous difference schemes. So after all I will start estimating the parameters values with a numerical data extracted, for instance, from a file. Thanks for your on-time help!

[Torsten]

The coeffcients c_0,...,c_n will be your unknown parameters.

Call lsqcurvefit with the vector of x_i as the xdata vector and the measured temperatures u_i as the ydata vector. Supply an initial guess of the vector of heat capacities c_i for lsqcurvefit.

lsqcurvefit will hand a suggestion for the parameter vector c_i of heat capacities to function "fun". Now in function fun, call the ODE/PDE solver of your choice to determine the temperatures u_i that result from the heat capacities c_i and return the vector of the calculated temperatures u_i to lsqcurvefit. Since the ODE/PDE solver will in general need heat capacities also in between the x_i, you will need to interpolate the c_i for a given x value in between x_0 and x_n (e.g. using Matlab's interp1).

If you already know a specific functional form for the heat capacity (e.g. c(u) =a + b*u + c*u^2)), the procedure will be easier. Now a,b and c will be the unknown parameter vector. Everything else will be the same as above, but you won't need to interpolate in between the x_i.

[Lyudmil Yovkov]

@Torsten, yes, I intend to search an approximation of the function c(u) as a quadratic one since some of the online resources about heat transfer recommend such a dependence. Thanks again!

[Lyudmil Yovkov]

I may have misunderstood you, this is why I'm attaching a source and a figure. It's something completely new for me and I need a bit of time to practice and get used to it.

%=================================
% u'(x) = f(x,u), 0 < x <= R, R = 4
% u(0) = u0 - initial condition
% u_{exact} = x^3 + 1 - exact solution
% f(x,u) = (3x^2) * (2x^6 + 4x^3 + 5) / (2u^2 + 3)
% We search c(u) = a + b*u + c*u^2.
%=================================
clear; clc;

% input data
xdata = linspace(0,4,200); 
ydata = xdata.^3 + 1;

% initial guess for p = [a,b,c]
p0 = [1,1,1];

% polynomial form of the
% approximate function c(u)
f = @(p,u) p(1) * u.^2 + ...
           p(2) * u + ...
           p(3);

% constraints for the parameters a, b, c
lb = [0.1,0,0];
ub = [3,3,3];

% calling lsqcurvefit
p = lsqcurvefit(f,p0,xdata,ydata,lb,ub);
c = @(u) p(1) * u.^2 + ...
         p(2) * u + ...
         p(3);
     
% right hand-side of the ODE
r_side = @(x,u) ... 
         3 * x.^2 * ...
         (2 * x.^6 + 4 * x.^3 + 5) ./ ...
         c(u);
% grid
x = xdata;
h = x(2) - x(1);
N = numel(x);

% exact solution
u = x.^3 + 1;

% approximate solution
y = zeros(1,N);
y(1) = u(1);

% Heun method
for i = 1 : N-1
    ym = y(i) + h * r_side(x(i),y(i));
    fl = r_side(x(i),y(i));
    fr = r_side(x(i+1),ym);
    y(i+1) = y(i) + (h/2) * (fl + fr);
end

% plot results
figure(1)
set(gcf,'color','w')
%%%
subplot(2,2,1)
plot(x,u,'b','LineWidth',3)
hold on, grid on
plot(x,y,'g--','LineWidth',3)
legend('\bf{Exact}','\bf{Approximated}')
set(gca,'FontSize',16)
%%%
subplot(2,2,2)
err = abs(y-u);
plot(x,err,'r:','LineWidth',3)
hold on, grid on
legend('\bf{|u_{exact}-u_{appr.}|}')
set(gca,'FontSize',16)
%%%
subplot(2,2,[3,4])
plot(u,c(u),'m','LineWidth',3)
hold on, grid on
legend('\bf{c(u)}')
title(['Estimated c(u) = ',num2str(p(1)),'u^2+',...
                 num2str(p(2)),'u+',...
                 num2str(p(3))])
set(gca,'FontSize',16)

[Torsten]

Your code tries to approximate x^3+1 in the interval [0;4] by a quadratic polynomial of the form p1*x^2+p2*x+p3. This is far off from what you want.

Proceed as follows:

Delete the definition of f in the upper part of your code and leave everything as is up to the call of lsqcurvefit.

In an extra .m-file, replace the definition of f by

function res = f(p,xdata)

fun = @(x,u) 3*x^2*(2*x^6+4*x^3+5)/(p(1)+p(2)*u+p(3)*u^2);

xspan = xdata;

u0 = 1;

[X,U] = ode15s(fun,xspan,u0);

res = U(:,1);

end

What is returned for p from lsqcurvefit ?

It should be approximately

p(1) = 3, p(2) = 0, p(3) = 2.

[Lyudmil Yovkov]

@Torsten, the function returns which is very, very close to the real values!

Yes, it's definitely working according to what I wanted to achieve! Thanks a lot, mate!

Kind regards!

[Torsten]

The result might even be better if you remove the setting of lower and upper bounds for the parameters.

[Lyudmil Yovkov]

@Torsten, is there any peculiarities if we wish to introduce dimensionless variables, for instance , where

[Torsten]

You mean

d/dx_cross = R * d/dx

?

[Lyudmil Yovkov]

Yes. I introduced the dimensionless variable and obtained what I wanted. But I've got a further question. Is there any difference in the method of estimation when we deal with an ODE and when we deal with a PDE?

[Torsten]

The structure of the program will be the same. In f, you will have to call the pdesolver instead of ode15s, and the computation time will be longer.

Answer 2

[Alan Stevens]

First plot the points, ui vs xi to see what sort of curve it might be. If it looks like it could be a polynomial (as in the case of u = x^2 + 1) then lookup help on polyfit.

[Lyudmil Yovkov]

@Alan Stevens, I gave quite a simple example in order to illustrate what I want to achieve. The real problem I've been dealing with is more complicated and nonlinear. In its formulation u(x) is a temperature and c(u) is a heat capacity. The heat capacity is usually a polynomial function of the temperature: For my goals it will be enough accurate to fit the function with a quadratic one. I've never before performed this inverse problem and this is why I decided to beg for help.

Answer 3

[Lyudmil Yovkov]

Hello to everyone,

I'm writing again for an advice. First let me introduce into what I've been working on. Given the Cauchy problem

where are known constants and is a known function of time. The functions have meaning of heat capacity and depend on the temperature T. For problem (1) an experimental data set as well is provided.

To solve (1) at known capacities and validate the results against the experimental data is a trivial procedure. I encounter difficulties when solving the inverse problem:

(2) Find the best fit for at known physical parameters and functions

I'm searching the capacity of the form , i.e. I'm trying to estimate the values of the best fit. There are some constraints imposed upon the values .

So thank to @Torsten, I successfully solved a modal problem of the type (2) but for test solutions prepared beforehand. Even at volatile constraints upon the values of the estimated parameters I obtain great results. I applied the ideas @Torsten suggested above to my real physics problem (2). It turned out that at small changes in the constraints or even at their full removal from the lsqcurvefit() procedure the results for also change significantly. Can you tell me why? Is the problem itself sensitive to the initial guess and the constraints, for instance?

I'm performing the experiments for pure water whose physical parameters are known. In particular, the function for pure water is known and I'm comparing my results with the analytical formula.

Best regards!

[Lyudmil Yovkov]

Can anybody advise me how to overcome the diffuculties when invoking lsqcurvefit() with different parameter constraints and explain to me the reason for the inaccurate numerical results? I'm attaching some plots.

Test 1, test 2, test 3 are performed on a test example with known function I'm comparing with. Here .

Fig. 1

Fig. 2

Fig. 3

Here are some of the results for the real physics problem I'm dealing with (Test 1 and test 2 below). Small changes, for instance, in the input constraints lead to significant displacement. The estimated -curve is not close to the exact one.

Fig. 4:

lb = [4.1, -2.3e-03, 4.8e-05, -4.8e-07, 1.7e-09],

ub = [4.3, -2.1e-03, 5.0e-05, -4.6e-07, 1.9e-09]

Fig. 5:

lb = [3.0, -3.0e-03, 3.5e-05, -5.0e-07, 1.0e-09],

ub = [5.0, -1.5e-03, 5.0e-05, -3.5e-07, 2.5e-09]

Thanks if you find a little time to respond to me and help me find the source of the errors.

Best regards!