Fitting the numerical solution of a PDE with one parameter to some data
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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.
0 Commenti
Risposta accettata
Torsten
il 22 Giu 2018
Modificato: Matt J
il 3 Mag 2023
The procedure can easily be adapted to work with a PDE instead of ODEs.
Best wishes
Torsten.
22 Commenti
Torsten
il 26 Giu 2018
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.
Più risposte (1)
Zana Taher
il 7 Apr 2019
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 Commenti
Torsten
il 20 Ago 2019
Modificato: Torsten
il 20 Ago 2019
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 ?
Vedere anche
Categorie
Scopri di più su Geometry and Mesh in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!