| 查看: 1262 | 回复: 5 | ||||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||||
[求助]
matlab lsqnonlin函数进行非线性动力学参数拟合，求高手指点. 已有1人参与
|
||||
|
function kinetics1 clear all clc a0=[10^11 81 0.3]; lb=[0 0 0];ub=[+inf +inf +inf]; c0=[0.0525564 100]; c16exp=... [0.0525564 100 0.0403533 73.3083 0.0267938 50.3759 0.015037 42.1053 0.0082474 24.812 0.00415801 13.9098 0.00278195 11.6541 0.00185881 10.1504 0.00183435 9.77444 0.0015879 9.77444]; c14exp=... [0.0525564 100 0.0493758 89.8496 0.0441658 81.203 0.0378303 72.9323 0.0357841 67.6692 0.0337337 61.6541 0.0276196 54.8872 0.0239832 47.7444 0.018766 40.6015 0.0160359 35.7143 0.0153455 34.5865]; c12exp=... [0.0526 100 0.0516 98.1203 0.0507 94.7368 0.0484 93.2331 0.0464 89.8496 0.0439 80.4511 0.0429 80.8271 0.0407 77.8195 0.0393 76.3158 0.0381 75.5639 0.0365 72.1805]; [a,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@objlm,a0,lb,ub,[],c0,c16exp,c14exp,c12exp); ci=nlparci(a,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为：\n'),a function f=objlm(a,c0,c16exp,c14exp,c12exp) %160度 tspan1=[0 0.1005 0.2143 0.3360 0.4365 0.5661 0.6667 0.7593 0.8677 0.9603]; [t,c16]=ode45(@kinetic16,tspan1,c0,[],a); f1=c16(:,1)-c16exp(:,1); f2=c16(:,2)-c16exp(:,2); %140度 tspan2=[0 0.1005 0.1984 0.2857 0.3571 0.4471 0.5529 0.6746 0.8042 0.9074 0.9683]; [t,c14]=ode45(@kinetic14,tspan2,c0,[],a); f3=c14(:,1)-c14exp(:,1); f4=c14(:,2)-c14exp(:,2); %120度 tspan3=[0 0.1005 0.2011 0.2937 0.4021 0.5000 0.5979 0.7090 0.8069 0.8862 0.9894]; [t,c12]=ode45(@kinetic12,tspan3,c0,[],a); f5=c12(:,1)-c12exp(:,1); f6=c12(:,2)-c12exp(:,2); f=[f1;f2;f3;f4;f5;f6]; %构造目标函数 function dc16dt=kinetic16(t,c16,a) %160度 T1=433.15;R=0.008314; dc16dt=... [-(a(1).*exp(-a(2)./(R.*T1)).*((c16(2)./100).^(a(3))).*c16(1).*0.943) -(7.*a(1).*exp(-a(2)./(R.*T1)).*((c16(2)./100).^(a(3))).*c16(1).*0.943)]; function dc14dt=kinetic14(t,c14,a) %140度 T2=413.15;R=0.008314; dc14dt=... [-(a(1).*exp(-a(2)./(R.*T2)).*((c14(2)./100).^(a(3))).*c14(1).*0.926) -(7.*a(1).*exp(-a(2)./(R.*T2)).*((c14(2)./100).^(a(3))).*c14(1).*0.926)]; function dc12dt=kinetic12(t,c12,a) %120度 T3=393.15;R=0.008314; dc12dt=... [-(a(1).*exp(-a(2)./(R.*T3)).*((c12(2)./100).^(a(3))).*c12(1).*0.907) -(7.*a(1).*exp(-a(2)./(R.*T3)).*((c12(2)./100).^(a(3))).*c12(1).*0.907)]; |
回复此楼» 猜你喜欢
综述论文作为代表作会不会影响评审专家的印象分？
已经有12人回复
-
2019年青年基金涵评意见，大家看看几个A，几个B?
已经有4人回复
-
应该是93bebmhtak前后十一个字符比较关键
已经有29人回复
-
奇怪，两个人的filecode固定段从头到尾一模一样
已经有9人回复
-
是这周出结果还是下周出结果？
已经有3人回复
-
小木虫上这么多卖论文的，真有人买论文么？感觉没必要啊
已经有9人回复
-
有时候，自然基金真的不能太认真 （我的申报经验）
已经有4人回复
-
帮忙看看fileCode
已经有13人回复
-
为什么网上很多人说本周 12号出结果
已经有7人回复
-
读博的好处
已经有3人回复
4楼2016-01-13 14:04:11
Mr__Right
专家顾问 (著名写手)
-
专家经验: +31 - 应助: 317 (大学生)
- 金币: 14456.3
- 散金: 500
- 红花: 54
- 帖子: 2716
- 在线: 950.6小时
- 虫号: 1972612
- 注册: 2012-09-04
- 性别: GG
- 专业: 应用数学方法
- 管辖: 程序语言
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
898766282: 金币+5, ★有帮助 2016-01-12 14:02:40
898766282: 金币+5, ★有帮助 2016-01-12 16:41:14
898766282: 金币+10, ★有帮助 2016-01-13 13:59:33
感谢参与，应助指数 +1
898766282: 金币+5, ★有帮助 2016-01-12 14:02:40
898766282: 金币+5, ★有帮助 2016-01-12 16:41:14
898766282: 金币+10, ★有帮助 2016-01-13 13:59:33
|
这类问题“尝试”的工作量大。只给一点提示： 1。拟合问题最终数值上都是“最优化”问题，这个是相对复杂的非线性非凸优化问题 2。这类优化问题的求解，常用的迭代方法或基于牛顿法延伸出来的方法，通常无法找到好的最优解 理论上，这类问题，只有唯一的全局最优解的情况下，没有数值算法能够确保必然收敛到这个点 现在通用的方法是： 先用全局优化算法（differential-evolution差分演化，最常用，粒子群PSO）作初步筛查，然后再用nonlincon，fmincon之类函数自带的局部优化算法 通常要多算几次做比较，可能还要调整优化的参数，直到找到最佳结果 |
2楼2016-01-12 13:03:33
3楼2016-01-12 16:41:04
Mr__Right
专家顾问 (著名写手)
-
专家经验: +31 - 应助: 317 (大学生)
- 金币: 14456.3
- 散金: 500
- 红花: 54
- 帖子: 2716
- 在线: 950.6小时
- 虫号: 1972612
- 注册: 2012-09-04
- 性别: GG
- 专业: 应用数学方法
- 管辖: 程序语言
5楼2016-01-13 19:12:34
