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

5 次查看（过去 30 天）

matnewbie 2018-6-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/406938-fitting-the-numerical-solution-of-a-pde-with-one-parameter-to-some-data

编辑： Matt J 2023-5-3

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.

采纳的回答

Torsten 2018-6-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/406938-fitting-the-numerical-solution-of-a-pde-with-one-parameter-to-some-data#answer_325814

编辑：Matt J 2023-5-3

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.

22 个评论
显示 20更早的评论隐藏 20更早的评论

matnewbie 2018-6-22

在 MATLAB Online 中打开

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 2018-6-26

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 2018-6-27

编辑：matnewbie 2018-7-2

Thanks for your precious help.

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

请先登录，再进行评论。

更多回答（1 个）

Zana Taher 2019-4-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/406938-fitting-the-numerical-solution-of-a-pde-with-one-parameter-to-some-data#answer_369474

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

% 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

4 个评论
显示 2更早的评论隐藏 2更早的评论

matnewbie 2019-8-20

That's not true, since the first experimental point (at

) is not reliable and you can discard it.

Torsten 2019-8-20

编辑：Torsten 2019-8-20

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 ?

请先登录，再进行评论。

请先登录，再回答此问题。