24小时热门版块排行榜    

查看: 1299  |  回复: 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的回帖

1135725495

铁杆木虫 (著名写手)

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

还有就是我按照你写的运行了,运行不了了。
出现下面的错误
Error: File: parafit.m Line: 72 Column: 1
Function definitions are not permitted in this context.
9楼2015-05-26 17:09:52
已阅   回复此楼   关注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的回帖
信息提示
请填处理意见