How can I estimate a coefficient of pde connecting pdepe and lsqcurvefit functions?

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.

Answer 1

[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