Fitting the numerical solution of a PDE with one parameter to some data

I need to fit the numerical solution of a PDE with one parameter to some data in MATLAB.

I already solved the PDE (by giving an arbitrary value to the parameter) and I have the data, but I am not sure if I can use nlinfit, since it requires an analytic function as input.

Another issue is how to pass the parameter to be fitted to pdepe to solve the PDE.

Answer 1

[Torsten]

https://de.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

The procedure can easily be adapted to work with a PDE instead of ODEs.

Best wishes

Torsten.

[matnewbie]

Thanks. It helped me, but I have still some issues... Here is the code:

DT0=1e-5;
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=...
lsqcurvefit(@MyPDE,DT0,radpos,radCprof);

function S=MyPDE(D_T,r,z)
m=1;
sol=pdepe(m,@pde,@ic,@bc,r,z);
% Extract the first solution component
S=sol(:,:,1);

      % Definition of the PDE
      function [c,f,s]=pde(r,z,c,DcDr)
          c=0.0152/D_T;
          f=DcDr;
          s=0;
      end

      % Initial condition
      function c0=ic(r)
          if r<0.5e-3
              c0=10;
          else
              c0=0;
          end
      end

      % Boundary conditions
      function [pl,ql,pr,qr]=bc(rl,cl,rr,cr,z)
          pl=0;
          ql=1;
          pr=0;
          qr=1;
      end
  end

The error I obtain is the following:

Not enough input arguments.

Error in TransverseDisp2>MyPDE (line 124) sol=pdepe(m,@pde,@ic,@bc,r,z);

Error in lsqcurvefit (line 202) initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in TransverseDisp2 (line 120) [B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@MyPDE,DT0,radpos,radCprof);

Caused by: Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

[Torsten]

"MyPde" expects 2 arguments, not 3.

[matnewbie]

I don't understand... Can you clarify by writing code?

[Torsten]

You write

function S=MyPDE(D_T,r,z)

with three arguments D_T,r,z.

The documentation of lsqcurvefit says that there have to be two arguments, not three.

[matnewbie]

Ok, but I have a PDE with a function (C) of two variables (r, z), not an ODE like in the example you mentioned.

So the extension to a PDE is not so trivial at the end...

[Torsten]

Convert the two vectors for r and z into a matrix (rz) using "meshgrid".

[matnewbie]

Ok, but in which point of the code? For sure not where the PDE is defined... Can you provide a working example?

[Torsten]

Before calling lsqcurvefit ; radpos and radCprof have to be (nr x nz) matrices, I guess. And S=sol(:,:,1) already is in this (nr x nz) format.

[matnewbie]

It doesn't work. I got the same error message, after I changed the code in this way:

r=radpos';
z=cellsize*(1:nx);

[rr,zz]=meshgrid(r,z);

DT0=1e-6;
[D,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@MyPDE,DT0,rr,zz);

[Torsten]

Generate an arbitrary (nz*nr)-matrix A.

Write your concentration measurements in a second (nzxnr)-matrix C.

Call lsqcurvefit as

[D,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@MyPDE,DT0,A,C);

Leave everything else as it is.

Best wishes

Torsten.

[matnewbie]

I tried your suggestion, but I got the same error message. Probably because now is the PDE that doesn't understand the two matrices.

This is the error message:

Not enough input arguments.

Error in TransverseDisp2>MyPDE (line 76) sol=pdepe(m,@pde,@ic,@bc,r,z);

Error in lsqcurvefit (line 202) initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in TransverseDisp2 (line 72) [D,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@MyPDE,DT0,Pmatrix,Cmatrix);

Caused by: Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

[Torsten]

DT0=1e-5;
A=...;
C=...;
r=...;
z=...;
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=...
lsqcurvefit(@(x,xdata)MyPDE(x,xdata,r,z),DT0,A,C);

function S=MyPDE(D_T,A,r,z)
m=1;
sol=pdepe(m,@pde,@ic,@bc,r,z);
% Extract the first solution component
S=sol(:,:,1);

        % Definition of the PDE
        function [c,f,s]=pde(r,z,c,DcDr)
            c=0.0152/D_T;
            f=DcDr;
            s=0;
        end
        % Initial condition
        function c0=ic(r)
            if r<0.5e-3
                c0=10;
            else
                c0=0;
            end
        end
        % Boundary conditions
        function [pl,ql,pr,qr]=bc(rl,cl,rr,cr,z)
            pl=0;
            ql=1;
            pr=0;
            qr=1;
        end
    end

[matnewbie]

Who are r and z in this case? Vectors?

[Torsten]

Yes, the usual spatial coordinate vectors.

[matnewbie]

Ok now it works, but least square fitting doesn't produce any valuable result. I got this output in the MATLAB console:

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point is less than the default value of the optimality tolerance.

stopping criteria details

[Torsten]

That's where the real work starts ...

Maybe varying the diffusion coefficient in a certain range, calculating the sum of differences squared between measured and simulated data for each case and choosing the coefficient with the smallest value for this sum could be worth a try.

[matnewbie]

I tried different random values (from 1e-3 to 1e-9) for the dispersion coefficient, but the message I got was always the same.

[Torsten]

"pdepe" must compute the concentration profiles very accurately for that lsqnonlin has a chance to calculate the derivatives of the concentrations with respect to the diffusion coefficient. That's why a very fine grid in r-direction and a small error tolerance for "pdepe" is essential.

[matnewbie]

And how could I refine the grid or put a small error tolerance for pdepe?

[Torsten]

Does S = sol(:,:,1) change significantly if you vary the diffusion coefficient ?

[Torsten]

And how could I refine the grid or put a small error tolerance for pdepe?

a) by using more points for the r-mesh

b) https://de.mathworks.com/help/matlab/ref/odeset.html (Variables "RelTol" and "AbsTol").

Best wishes

Torsten.

[matnewbie]

Thanks for your precious help.

In any case I noticed that also nlinfit works, so I will use it instead of lsqcurvefit.

Answer 2

[Zana Taher]

Using the matlab function lsqnonlin will work better than lsqcurvefit. Below is an example code for using lsqnonlin with pdepe to solev parameters for richard's equation.

global exper

exper=[-0.4,-112,-121,-128,-131,-131.1,-133.2,-133.1,-134.4,-134.7,-135.12,-135.7,-135.8,-135.1,-135.4,-135.8,-135.9,-134,-135.7,-135.1,-135.3,-135.4,-136,-139,-136];

% p=[0.218,0.52,0.0115,2.03,31.6];

% qr f a n ks

p0=[0.218,0.52,0.0115,2.03,31.6];

lb=[0.01,0.4,0,1,15]; %lower bound

ub=[Inf,Inf,10,10,100]; %upper bound

options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter');

[p,resnorm,residual,exitflag,output] = lsqnonlin(@richards1,p0,lb,ub,options)

uu=richards1(p);

global t

plot(t,uu+exper')

hold on

plot(t,exper,'ko')

function uu=richards1(p)

% Solution of the Richards equation

% using MATLAB pdepe

%

% $Ekkehard Holzbecher $Date: 2006/07/13 $

%

% Soil data for Guelph Loam (Hornberger and Wiberg, 2005)

% m-file based partially on a first version by Janek Greskowiak

%--------------------------------------------------------------------------

L = 200; % length [L]

s1 = 0.5; % infiltration velocity [L/T]

s2 = 0; % bottom suction head [L]

T = 4; % maximum time [T]

% qr = 0.218; % residual water content

% f = 0.52; % porosity

% a = 0.0115; % van Genuchten parameter [1/L]

% n = 2.03; % van Genuchten parameter

% ks = 31.6; % saturated conductivity [L/T]

x = linspace(0,L,100);

global t

t = linspace(0,T,25);

options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7)

u = pdepe(0,@unsatpde,@unsatic,@unsatbc,x,t,options,s1,s2,p);

global exper

uu=u(:,1,1)-exper';

% figure;

% title('Richards Equation Numerical Solution, computed with 100 mesh points');

% subplot (1,3,1);

% plot (x,u(1:length(t),:));

% xlabel('Depth [L]');

% ylabel('Pressure Head [L]');

% subplot (1,3,2);

% plot (x,u(1:length(t),:)-(x'*ones(1,length(t)))');

% xlabel('Depth [L]');

% ylabel('Hydraulic Head [L]');

% for j=1:length(t)

% for i=1:length(x)

% [q(j,i),k(j,i),c(j,i)]=sedprop(u(j,i),p);

% end

% end

%

% subplot (1,3,3);

% plot (x,q(1:length(t),:)*100)

% xlabel('Depth [L]');

% ylabel('Water Content [%]');

% -------------------------------------------------------------------------

function [c,f,s] = unsatpde(x,t,u,DuDx,s1,s2,p)

[q,k,c] = sedprop(u,p);

f = k.*DuDx-k;

s = 0;

end

% -------------------------------------------------------------------------

function u0 = unsatic(x,s1,s2,p)

u0 = -200+x;

if x < 10 u0 = -0.5; end

end

% -------------------------------------------------------------------------

function [pl,ql,pr,qr] = unsatbc(xl,ul,xr,ur,t,s1,s2,p)

pl = s1;

ql = 1;

pr = ur(1)-s2;

qr = 0;

end

%------------------- soil hydraulic properties ----------------------------

function [q,k,c] = sedprop(u,p)

qr = p(1); % residual water content

f = p(2); % porosity

a = p(3); % van Genuchten parameter [1/L]

n = p(4); % van Genuchten parameter

ks = p(5); % saturated conductivity [L/T]

m = 1-1/n;

if u >= 0

c=1e-20;

k=ks;

q=f;

else

q=qr+(f-qr)*(1+(-a*u)^n)^-m;

c=((f-qr)*n*m*a*(-a*u)^(n-1))/((1+(-a*u)^n)^(m+1))+1.e-20;

k=ks*((q-qr)/(f-qr))^0.5*(1-(1-((q-qr)/(f-qr))^(1/m))^m)^2;

end

end

end

[matnewbie]

Thank you for your suggestion. How could I adapt my PDE to your solution strategy?

The PDE with the boundary conditions is the following:

I tried to use the following form for the PDE:

function S=MyPDE(D_T,r,z)
global Uavg
m=1;
sol=pdepe(m,@pde,@ic,@bc,r,z);
% Extract the first solution component
S=sol(:,:,1);
      % Definition of the PDE
      function [c,f,s]=pde(r,z,c,DcDr)
          c=Uavg/D_T;
          f=DcDr;
          s=0;
      end
      % Initial condition
      function c0=ic(r)
          if r<0.5e-3
              c0=6;
          else
              c0=0;
          end
      end
      % Boundary conditions
      function [pl,ql,pr,qr]=bc(rl,cl,rr,cr,z)
          pl=0;
          ql=1;
          pr=0;
          qr=1;
      end
end

but I obtain a mismatch for the concentration values, even if the transverse Peclet number (related to D_T, which I put equal to a value of 1e-5 without performing lsqnonlin, just to have an idea of the solution) has a reasonable value of 12, as you can see from the following plots:

[Torsten]

You'll never get a profile from the PDE solution with a maximum concentration in the interior of your domain. Thus your model equation seems to be inadequate for the experiment you performed.

[matnewbie]

That's not true, since the first experimental point (at ) is not reliable and you can discard it.

[Torsten]

You set dc/dr = 0 at both ends of your domain. So no mass leaves the system. Thus the initial mass of your species can be obtained by integrating c from r=0 to r=R. If you do this for both experiment and model, it is obvious that the initial mass of the species in the system for the simulation must have been much lower than in your experiment. So how can you expect that the curves agree ?