Problems in passing data from lsqcurvefit to pdepe solver
2 次查看(过去 30 天)
显示 更早的评论
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.
回答(1 个)
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
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!