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

5 次查看（过去 30 天）
matnewbie
matnewbie 2018-6-22
编辑： 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
Torsten 2018-6-22
编辑：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更早的评论
Torsten
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
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
Zana Taher 2019-4-7
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
  4 个评论
显示 2更早的评论
matnewbie
matnewbie 2019-8-20
That's not true, since the first experimental point (at ) is not reliable and you can discard it.
Torsten
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 ?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by