| 查看: 134 | 回复: 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 ] |
回复此楼» 猜你喜欢
有多少人是今天查系统知道结果的？
已经有19人回复
-
国社科又开始会评了，不知道这次命运如何
已经有9人回复
-
能否申诉？
已经有5人回复
-
为什么到现在没收到通知？
已经有5人回复
-
有没有仍没收到信息的
已经有8人回复
-
面上意见出来了
已经有18人回复
-
找导师
已经有7人回复
-
29号明天会评吗
已经有3人回复
-
麻烦专家们看看评委们的意见（F口面上）
已经有11人回复
-
我就是申请一个面上项目而已，这评审意见是按照杰青的条件评的吧？
已经有11人回复
投票: