| 查看: 130 | 回复: 0 | |||
| 当前主题已经存档。 | |||
| 【有奖交流】积极回复本帖子，参与交流，就有机会分得作者 Shoney 的 30 个金币 | |||
Shoney木虫 (著名写手)
|
[交流]
【求助】大家帮忙看看这个求最优参数的程序问题在哪儿？
|
||
|
是一个根据实验i数据求最优扩散系数的问题，即有一种物质扩散入一个圆柱体内，实验测得各个时间点的浓度，方程是一个PDE方程。运行结果没提示错误，但结果就是不对，最优的Deff值老是等于初值，不知道问题出在哪儿。（附件是源程序）程序如下： function ACdisc_Diffusion clear all;clc; global Cm; %动力学数据 tlist=[0 1 2 5 10 20 30 60]; %测量的时间点，单位min Cm=[0 4.30 5.20 6.17 7.53 8.41 8.80 9.04]; %测得的AC中DMSO平均浓度，单位%w/w %非线性拟合 C0=0; Deff0=0.0116; %单位mm2/min [Deff,resnorm,resid,exitflag,output,lambda,jacobian]=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; %方程的系数 %网格化（三角形网格划分及网格细化） [p,e,t]=initmesh(g); [p,e,t]=refinemesh(g,p,e,t); [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: [1x137 char] lambda = lower: 0 upper: 0 jacobian = All zero sparse: 8-by-1 ci = 1.0e+007 * -9.8973 9.8973 [ Last edited by Shoney on 2009-9-15 at 10:42 ] |
回复此楼» 猜你喜欢
好奇怪的filecode
已经有5人回复
-
应该是93bebmhtak前后十一个字符比较关键
已经有21人回复
-
综述论文作为代表作会不会影响评审专家的印象分？
已经有10人回复
-
奇怪，两个人的filecode固定段从头到尾一模一样
已经有3人回复
-
国自然结果
已经有5人回复
-
8月时间戳变的，举个手。玩一下，释放压力
已经有12人回复
-
我的国基提前知道中了，可是同事的操作让我实在接受不了，怎么会有这样的人
已经有8人回复
-
据悉今年马上要出结果了
已经有8人回复
-
关于代码变化问题，想知道的进来
已经有23人回复
-
FileCode能看出啥？
已经有6人回复
20