24小时热门版块排行榜    

查看: 1298  |  回复: 10
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

1135725495

铁杆木虫 (著名写手)

[求助] 求大家帮我看看下面中matlab中的程序问题出在哪里?已有1人参与

求动力学参数
方程为:r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
t          c
0        8.769713879
3        5.909461709
6        5.736047766
9        4.756936544
12        4.505799773
15        4.354160734
18        4.124025932
21        3.994725485
31        3.732949397
41        3.680613046
51        3.537445782
61        3.283577643
求k1,k2,k3,,k4
下面是写的程序,大家帮忙看看,为什么结果是直线呢?
function parafit
%  
% r1 =(k(1)*k(2)*c^2-0.25*k(3)*k(4)*(9.404-c)^2)/(1+k(2)*c+0.5*k(4)*(9.404-c));
%
% dCAdt = - r1;
clear all
clc
%        t/min   CA     / mol/L
  Kinetics=[0        7.029559717
1        3.960706476
2        3.163862085
3        2.963057756
6        2.820434159
9        2.698109429
12        2.556212898
15        2.447100967
18        2.207246045
21        2.13921916
31        2.000403189
41        1.917172016
51        1.851299168
61        1.636287388];
c=Kinetics( :,1)';
k0 = [0.671   0.02253  0.641   0.04028];         % 参数初值
lb = [0 0 0 0];                   % 参数下限
ub = [100 100  100  100];    % 参数上限
x0 = [7.029559717];
yexp = Kinetics;                  
warning off
% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
warning off
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
warning off
% 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('ca-pr','ca-real')
x(:,1)
figure;plot(t,x);
hold on

% ------------------------------------------------------------------
function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
f =  sum((y(:,2)-yexp(:,2)).^2) ;
% ------------------------------------------------------------------
function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0 1 2 3 6 9 12 15 18 21 31 41 51 61];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);

f1 = y(:,2) - yexp(:,2);

f = [f1];
% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
Kinetics=[0        7.029559717
1        3.960706476
2        3.163862085
3        2.963057756
6        2.820434159
9        2.698109429
12        2.556212898
15        2.447100967
18        2.207246045
21        2.13921916
31        2.000403189
41        1.917172016
51        1.851299168
61        1.636287388];
c=Kinetics(:,1)';
for i=1:14
    c(i)=c(:,i);
   
r1 =(k(1)*k(2)*c(i)^2-0.25*k(3)*k(4)*(9.404-c(i))^2)/(1+k(2)*c(i)+0.5*k(4)*(9.404-c(i)));
dCAdt = - r1;
dxdt = [dCAdt];
end
结果如下:

求大家帮我看看下面中matlab中的程序问题出在哪里?
H8_DJFZ~6NZC}07J0EXL1]L.png
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)

引用回帖:
7楼: Originally posted by 1135725495 at 2015-05-26 17:05:15
谢谢啊!请问这是改好的吗?我可以直接用吗?...

改好,但是初值的选取不当,拟合结果不是很好。你可以试试

» 本帖已获得的红花(最新10朵)

The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
8楼2015-05-26 17:08:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 11 个回答

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
1135725495: 金币+10, 有帮助 2015-05-26 15:35:53
两组参考结果:

Correlation Coef. (R): 0.986609345793065
R-Square: 0.97339800120622
Determination Coef. (DC): 0.973395608126228
F-Statistic: 73.3952829093829
Parameters Best Estimate
-------------------- -------------
k1 0.0633244071641639
k2 -0.123667552854961
k3 0.0140433904054541
k4 -0.681989563148558
====== Output Results =====
File: Data file - 1
No Obs. c Cal. c
1 5.909461709 6.03887691479301
2 5.736047766 5.36051185683508
3 4.756936544 4.91507500618295
4 4.505799773 4.59634816701972
5 4.354160734 4.35722735701414
6 4.124025932 4.17236164824689
7 3.994725485 4.02646466838814
8 3.732949397 3.71546906582311
9 3.680613046 3.55315135651946
10 3.537445782 3.46435427377623
11 3.283577643 3.41460522937653

Correlation Coef. (R): 0.98660934579373
R-Square: 0.973398001207533
Determination Coef. (DC): 0.973395608126228
F-Statistic: 73.3952829093831
Parameters Best Estimate
-------------------- -------------
k1 0.481614078106072
k2 1.61663483596038
k3 0.0205007668163867
k4 46.4476555663092
====== Output Results =====
File: Data file - 1
No Obs. c Cal. c
1 5.909461709 6.03887691237705
2 5.736047766 5.36051185622609
3 4.756936544 4.91507500636805
4 4.505799773 4.5963481675709
5 4.354160734 4.35722735771124
6 4.124025932 4.17236164896553
7 3.994725485 4.02646466905497
8 3.732949397 3.71546906608907
9 3.680613046 3.55315135632873
10 3.537445782 3.46435427319561
11 3.283577643 3.41460522849613
MATLAB、MS小问题、普通问题请发帖求助!时间精力有限,恕不接受无偿私信求助。
2楼2015-05-26 13:56:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

1135725495

铁杆木虫 (著名写手)

引用回帖:
2楼: Originally posted by 月只蓝 at 2015-05-26 13:56:42
两组参考结果:

Correlation Coef. (R): 0.986609345793065
R-Square: 0.97339800120622
Determination Coef. (DC): 0.973395608126228
F-Statistic: 73.3952829093829
Parameters Best Estimate
--------- ...

您好,您的计算结果我看到了,我想说的是我上面编的Matlab程序正确吗?为什么会出现那样的计算结果?还有就是您计算的方法是最小二乘吗?我现在需要用matlab计算,因为写论文需要。所以您能不能帮忙看一下程序是否正确?
3楼2015-05-26 15:38:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

引用回帖:
3楼: Originally posted by 1135725495 at 2015-05-26 15:38:21
您好,您的计算结果我看到了,我想说的是我上面编的Matlab程序正确吗?为什么会出现那样的计算结果?还有就是您计算的方法是最小二乘吗?我现在需要用matlab计算,因为写论文需要。所以您能不能帮忙看一下程 ...

拟合曲线和散点差异明显,可能是因为参数初值选取得不好。
MATLAB、MS小问题、普通问题请发帖求助!时间精力有限,恕不接受无偿私信求助。
4楼2015-05-26 16:25:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
信息提示
请填处理意见