大家帮忙看看这个求最优参数的程序问题在哪儿？
本帖最后由 shoney 于 2009-9-15 10:40 编辑是一个根据实验i数据求最优扩散系数的问题，即有一种物质扩散入一个圆柱体内，实验测得各个时间点的浓度，方程是一个PDE方程。运行结果没提示错误，但结果就是不对，最优的Deff值老是等于初值，不知道问题出在哪儿。（附件是源程序）。程序如下：
function ACdisc_Diffusion
clear all;clc;
global Cm;
%动力学数据
tlist=; %测量的时间点，单位min
Cm=; %测得的AC中DMSO平均浓度，单位%w/w
%非线性拟合
C0=0;
Deff0=0.0116; %单位mm2/min
=lsqnonlin(@OptObjFunc,Deff0,[],[],[],tlist,C0,Cm)
ci=nlparci(Deff,resid,jacobian)
%--------------------------------------------------------------------------
%求解PDE
function C=OptObjFunc(Deff,tlist,C0,Cm)
%PDE问题定义及参数初始化
g='ACdiscg1'; %定义求解域，几何尺寸单位mm
b='ACdiscb1'; %定义边界条件
d=['x/' num2str(Deff)];c='x';a=0;f=0; %方程的系数
%网格化（三角形网格划分及网格细化）
=initmesh(g);
=refinemesh(g,p,e,t);
=refinemesh(g,p,e,t);
p=jigglemesh(p,e,t);
%绘制PDE三角形网格
pdemesh(p,e,t);
u=parabolic(C0,tlist,b,p,e,t,c,a,f,d); %计算得到的AC中任意时间任意位置DMSO浓度
%积分求解DMSO平均浓度的计算值
xi=linspace(0,3,100);yi=linspace(0,0.38,100);
Cc2(1,1)=0; %赋值0，因进行网格转化时，值由线性插值得到；结果Cc2(1,1)≠0(边界值=9.04)
for i=2:length(tlist)
Cc{1,i}=tri2grid(p,t,u(:,i),xi,yi); %各个时间点分别进行网格转化
Cc0=Cc{1,i};
Cc1(1,i)=dblquad(@Func,0,3,0,0.38,[],@quadl,xi,yi,Cc0); %积分求各个时间点对应的总浓度
Cc2(1,i)=Cc1(1,i)/(3*0.38); %平均浓度
end
C=Cc2-Cm; %残差
%--------------------------------------------------------------------------
function f=Func(x,y,xi,yi,Cc0)
f=interp2(xi,yi,Cc0,x,y);
运行结果如下：
Optimization terminated: first-order optimality less than OPTIONS.TolFun,
and no negative/zero curvature detected in trust region model.
Deff =
0.0116
resnorm =
2.7230
resid =
0 -1.0917 -0.7775 0.4317 0.6481 0.5194 0.2249 0.0000
exitflag =
1
output =
firstorderopt: 0
iterations: 0
funcCount: 2
cgiterations: 0
algorithm: 'large-scale: trust-region reflective Newton'
message:
lambda =
lower: 0
upper: 0
jacobian =
All zero sparse: 8-by-1
ci =
1.0e+007 *
-9.8973 9.8973
觉得这句d=['x/' num2str(Deff)];有问题。另外，你这个问题最好分开算，先解方程，然后再拟合。 觉得这句d=['x/' num2str(Deff)];有问题。另外，你这个问题最好分开算，先解方程，然后再拟合。
messenger 发表于 2009-9-15 11:26 http://forum.simwe.com/images/common/back.gif
一直就很讨厌MATLAB的这种符号计算的表达方式，很别扭、生涩。
建议用comsol来做你的问题，matlab自带的那个pde tool差太多了。 觉得这句d=['x/' num2str(Deff)];有问题。另外，你这个问题最好分开算，先解方程，然后再拟合。
messenger 发表于 2009-9-15 11:26 http://forum.simwe.com/images/common/back.gif这句没有问题，确定。因为单独解PDE时（Deff用估计值，含这句）没有问题 觉得这句d=['x/' num2str(Deff)];有问题。另外，你这个问题最好分开算，先解方程，然后再拟合。
messenger 发表于 2009-9-15 11:26 http://forum.simwe.com/images/common/back.gif按照您“分开算”的方法，我试了一下，自己一个一个的带入Deff值，然后比较残差平方和，在Deff=0.0122时，得到最小的残差平方和。但是当我把这个Deff值打入这个程序的时候，结果还是不对，郁闷啊……自己感觉lsqnonlin()这个函数也没用错啊 本帖最后由 huchengyao7890 于 2009-9-16 12:35 编辑
function ACdisc_Diffusion2
clear all;clc;
global y1;
t2=;
y1=;
beta0=;
=lsqnonlin(@fun446,beta0,[],[],[],t2,y1)
% --------------------------------------------------------------
function F=fun446(beta,t2,y1)
m = 0;
x = ;
t = ;
sol = pdepe(m,@pdex5pde,@pdex5ic,@pdex5bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
nt = length(t2);
iok = 2:nt;
for j = iok
= pdeval(m,x,u2(j,:),0);
end
F=I+beta(1)+beta(2).*t2-y1;
% --------------------------------------------------------------
function = pdex5pde(x,t,u,DuDx)
c = ;
f = .* DuDx;
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F];
% --------------------------------------------------------------
function u0 = pdex5ic(x);
u0 = ;
% --------------------------------------------------------------
function = pdex5bc(xl,ul,xr,ur,t)
pl = ;
ql = ;
pr = ;
qr = ;
谁帮我看看我的问题啊，我是直接用pdepe求解的一维偏微分方程的，可是好像无法计算，帮帮我吧，一运行总说Error using ==> beta
Not enough input arguments.，可是我两个参数，两个初始值，没有问题啊。 我向老师请教过，我估计是你的参数根本就没有传递给你的偏微分方程，你只是写到系数里，可是求解偏微分方程的程序里并没有传递这个参数进去，不知道我说的对不对啊？ 我向老师请教过，我估计是你的参数根本就没有传递给你的偏微分方程，你只是写到系数里，可是求解偏微分方程的程序里并没有传递这个参数进去，不知道我说的对不对啊？
huchengyao7890 发表于 2009-9-16 15:34 http://forum.simwe.com/images/common/back.gif"可是求解偏微分方程的程序里并没有传递这个参数进去"这个怎么理解呢？Deff在PDE方程的系数里，那计算PDE的时候肯定是包含进了Deff啊？
页:
[1]