Problems in passing data from lsqcurvefit to pdepe solver

 Hi

I'm trying to fit a function to a three dimensional dataset of concentration of drug through a tumor by using lsqcurvefit.

The diffusion of the drug throught he tumor is being modeled by the partial differential equation

dC/dt=D/r^2*d/dr(r^2*dC/dr)-ke*C

D=diffusion r=radius of tumor ke=rate of elimination of drug C=concentration

And the parameters I am looking to fit are my D, ke, and co where C(r,t=0)=co

My boundary conditions are

dC/dr(r=0,t)=0 C(r=R,t)=Cr=1

I've written a pdepe solver function that recieves the data and my initial parameter guesses and another program which calls this function using lsqcurvefit. However, when I run the code I get the following errors for the guess of

D=8.64e-2; ke=8.64e-2; co=1;

Warning: Failure at t=3.482418e-002.  Unable to meet
integration tolerances without reducing the step
size below the smallest value allowed
(1.110223e-016) at time t. 
> In ode15s at 753
  In pdepe at 320
  In pset7p5_449 at 22
  In lsqcurvefit>objective at 266
  In optim\private\snls at 352
  In optim\private\lsqncommon at 153
  In lsqcurvefit at 258
  In pset7p5_449fit at 13
Warning: Time integration has failed. Solution is
available at requested time points up to
t=0.000000e+000. 
> In pdepe at 326
  In pset7p5_449 at 22
  In lsqcurvefit>objective at 266
  In optim\private\snls at 352
  In optim\private\lsqncommon at 153
  In lsqcurvefit at 258
  In pset7p5_449fit at 13
??? Error using ==> minus
Matrix dimensions must agree.

Error in ==> lsqcurvefit>objective at 267
         F = F - YDATA;

Error in ==> snls at 352
            newfvec =
            feval(funfcn{3},xcurr,varargin{:});

Error in ==> lsqncommon at 153
    [xC,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msg]=...

Error in ==> lsqcurvefit at 258
[xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB]
= ...

Error in ==> pset7p5_449fit at 13
[b resnorm residual
exitflag]=lsqcurvefit(@pset7p5_449,b0,data,C(:,1:8));

Oddly, for some guesses ([0 0 1] the code seems to be able to go on long enough to produce an estimate of my parameters but terminates early due to either the number of iterations or the TolX, which I have modulated to try to improve my results to no avail. I don't understand why the code is so sensitive to the initial guesses. It seems like the true parameters are likely to be [0.008 0.035 .0077].

My code is as follows

clear
load tumor_data.mat

beta=[.0432 .864 1];
otherguess=[0 0 1];
Dg=8.64e-2; kg=8.64e-2; c0=1;
b0=[Dg kg c0];
data.r=r;
data.t=t./86400;
options=optimset('MaxFunEvals',5000,'TolX',1e-6,'TolFun',1e-10);

[b resnorm residual exitflag]=lsqcurvefit(@pset7p5_449,b0,data,C);

Csol=pset7p5_449(b,data);

figure(1)
mesh(data.r,data.t,Csol)
xlabel('Distance in cm')
ylabel('Time in Days')
zlabel('Concentration')
%title(strcat('D=',num2str(D),' ke=',num2str(ke)))
axis tight

--------------------------------------------

function sol=pset7p5_449(b,data)

% xf = 1;
% nx = 200;
% tf = 7; %days
% nt = 200;

xmesh=data.r;
tspan=data.t;
options=[];
m=2;

% ke=8.64e-1; %1/day
% D=4.32e-2; %cm^2/day
Cr=1;

D=b(1);
ke=b(2);
Co=b(3);

sol=pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan,options);

function [c,f,s]=pdefun(x,t,u,DuDx)
    c = 1;
    f = D*DuDx;
    s = -ke*u;
end

function u0 = icfun(x)
    u0 = Co;
end

function [pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
    pl=0;
    ql=1;
    pr=ur-Cr;
    qr=0;
end

end

---------------

the data in the tumor fileset is for r, t, and C where r is a 1x20 t is a 1x8 and c is a 8x20

ANY help on this would be appreciated. I've been fidgeting with this to no avail for some time now.

[Walter Roberson]

Can you attach the .mat so people can test ?

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