How can I estimate a coefficient of pde connecting pdepe and lsqcurvefit functions?
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
I have similar problem to this one: http://www.mathworks.com/matlabcentral/answers/5603-problems-in-passing-data-from-lsqcurvefit-to-pdepe-solver
I am trying to fit only D coefficient and I am getting a following error:
??? Attempted to access fL(1); index out of bounds because numel(fL)=0.
Error in ==> pdepe>pdeodes at 363 up(j,1) = sL(j) + (m+1) * fL(j) / xi(1);
Error in ==> odearguments at 98 f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
I would be grateful for any help.
0 comentarios
Respuestas (1)
Zana Taher
el 7 de Abr. de 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
0 comentarios
Ver también
Categorías
Más información sobre Visualization en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!