| 查看: 1233 | 回复: 4 | |||
[求助]
求助matlab拟合动力学K值 已有1人参与
|
|
动力学模型为 dC(1)/dt= -((k(1)+k(2)))*C(1) dC(2)/dt= k(1)*C(1)-k(3)*C(2)-k(4)*C(2) dC(3)/dt= k(2)*C(1)-k(5)*C(3) dC(4)/dt= -k(6)*C(4)+k(3)*C(2) dC(5)/dt= k(4)*C(2)+k(5)*C(3)-k(7)*C(5) dC(6)/dt= k(6)*C(4)+k(7)*C(5) ; 实验数据 % t C(1) C(2) C(3) C(4) C(5) C(6) % ExpData =[ 0 238.5300 0 0 0 0 0 22.0000 77.1500 8.3600 13.9900 0 0.5100 0 36.0000 58.4500 13.0200 27.5800 0.0300 0.9200 0 55.0000 33.0900 16.8200 48.3100 0.0500 1.7100 0.0100 76.0000 15.9600 25.5900 55.2300 0.1200 3.0500 0.0500 96.0000 1.5500 25.7900 63.3900 0.4200 8.8400 0 115.0000 0.0400 19.1900 46.9500 1.2800 32.4400 0.1000 134.0000 0 17.8700 25.0700 2.0000 54.8500 0.2100 154.0000 0.0200 6.2800 15.7100 2.5300 74.9700 0.4900 173.0000 0 1.6900 3.0300 2.9600 91.3100 0.9900 177.0000 0 1.1000 1.3300 2.9300 93.6800 0.9500 178.0000 0 0.5400 0.3600 3.1000 94.9000 1.0800 188.0000 0 0.3700 0.5000 3.0400 94.7400 1.3500 ] 哪位大神能帮忙拟合一下k值呢？为什么我一直得不出结果  |
回复此楼» 猜你喜欢
B区学生调剂-兰州交通大学材料科学与工程学院
已经有12人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有13人回复
-
物理学I论文润色/翻译怎么收费?
已经有78人回复
-
桂林理工大学物理学专业招收调剂
已经有18人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
-
津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生一名
已经有0人回复
-
【原创讨论】从电子约束到物质编辑：一套可迭代的环形磁场科技树
已经有0人回复
-
【方案分享】单环磁场+轴心控制+偏转导出电子束约束系统（可行性实验）
已经有6人回复
-
【修正版】单环用磁约束低速电子实验方案（简化版）
已经有0人回复
-
桂林理工大学物理学专业招收调剂，还有三个名额！！！
已经有12人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab拟合曲面方程
已经有5人回复
- 使用MATLAB拟合动力学参数问题
已经有18人回复
- 用MATLAB拟合如下数据（求程序代码和拟合结果）
已经有12人回复
- matlab 非线性拟合
已经有16人回复
- Matlab同时拟合活化能，反应级数，速率常数
已经有26人回复
- 使用matlab最优化方法拟合获得多个动力学参数中的问题
已经有4人回复
- 求助用matlab拟合动力学方程
已经有9人回复
- 求助用matlab拟合动力学方程
已经有13人回复
- matlab 拟合求助
已经有8人回复
- 紧急求助，利用Matlab对实验数据进行拟合求解参数。
已经有27人回复
- matlab数据拟合问题求助
已经有7人回复
- matlab 拟合反应动力学参数结果很差。大家帮忙看一下
已经有14人回复
- matlab如何编写共享参数拟合程序
已经有8人回复
- 【求助】用matlab最优化方法进行参数拟合
已经有17人回复
- 【求助】拟合动力学方程求助
已经有13人回复
- 【求助】催化反应动力学matlab计算各基元反应的速率常数时，该如何避免较小量被忽略？
已经有3人回复
- 【求助】matlab曲面拟合表达式
已经有10人回复
- 【求助】使用MATLAB怎么实现拟合力场参数的程序？【已解决】
已经有15人回复
- 【求助】使用Matlab预估动力学方程问题
已经有13人回复
【答案】应助回帖
感谢参与，应助指数 +1
|
你这个是怎么思路做的？我这里有一个简单的程序，你可以看一眼。 function k1k2k3 format long clear all clc tspan = [0 30 50 80 140 170 200 230 260 290 320 360 400 460 520 580 600]; x0 = [9; 0]; k0 = [0.001 0 0.08]; lb = [0 0 0]; ub = [1 1 1]; data=... [ 30 7.939 1.458 50 7.687 1.535 80 7.289 1.602 140 6.658 1.717 170 6.531 1.722 200 6.218 1.671 230 5.979 1.620 260 5.591 1.550 290 5.414 1.488 320 4.968 1.433 360 4.692 1.350 400 4.438 1.319 460 4.144 1.294 520 4.041 1.294 580 4.052 1.287 600 4.030 1.292 ]; yexp = data(:,2:3); [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3)) fprintf(' The sum of the squares is: %.9e\n\n',resnorm) function f = ObjFunc(k,tspan,x0,yexp) % 目标函数 [t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k); Xsim1=Xsim(:,1); Xsim2=Xsim(:,2); ysim(:,1) = Xsim1(2:end); ysim(:,2) = Xsim2(2:end); %size(ysim(:,1)); %size(ysim(:,2)); %size(yexp(:,1)); %size(yexp(:,2)); f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2))]; function dCdt = KineticsEqs(t,C,k) % ODE模型方程 dCAdt = -k(1)*C(1)+k(2)*C(2); dCBdt = k(1)*C(1)-k(2)*C(2)-k(3)*C(2); dCdt = [dCAdt; dCBdt]; |
2楼2015-06-26 19:41:28
|
本帖内容被屏蔽 |
3楼2015-06-27 12:34:16
|
这个是我照书上程序改的，可是计算不出结果。 function KineticsEst5 clear all clc k0 = [0.5 0.5 0.5 0.5 0.5 0.5 0.5]; % 参数初值 lb = [0 0 0 0 0 0 0]; % 参数下限 ub = [+inf +inf +inf +inf +inf +inf +inf]; % 参数上限 C0 = [ 238.53 0 0 0 0 0 ]; ExpData =[0 238.53 0 0 0 0 0 22 77.15 8.36 13.99 0 0.51 0 36 58.45 13.02 27.58 0.03 0.92 0 55 33.09 16.82 48.31 0.05 1.71 0.01 76 15.96 25.59 55.23 0.12 3.05 0.05 96 1.55 25.79 63.39 0.42 8.84 0 115 0.04 19.19 46.95 1.28 32.44 0.1 134 0 17.87 25.07 2 54.85 0.21 ] t= ExpData(:,1)' yexp = ExpData(:,2:7); % yexp: 实验数据[x1 x4 x5 x6] % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],C0,yexp); fprintf('\n使用函数fmincon()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf('\tk3 = %.4f\n',k(3)) fprintf('\tk4 = %.4f\n',k(4)) fprintf('\tk5 = %.4f\n',k(5)) fprintf('\tk5 = %.4f\n',k(6)) fprintf('\tk5 = %.4f\n',k(7)) fprintf(' The sum of the squares is: %.1e\n\n',fval) k_fmincon = k; % 使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],C0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') Output % 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0 = k_fmincon; [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],C0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n') Output % ------------------------------------------------------------------ function f = ObjFunc4Fmincon(k,C0,yexp) tspan = [0 22 36 55 76 96 115 134 ]; [t C] = ode45(@KineticEqs,tspan,C0,[],k); y(:,1) = C(:,1); y(:,2:6) = C(:,2:6); f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2) ... + sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2) ... + sum((y(:,5)-yexp(:,5)).^2) + sum((y(:,6)-yexp(:,6)).^2); % ------------------------------------------------------------------ function f = ObjFunc4LNL(k,C0,yexp) tspan = [0 22 36 55 76 96 115 134 ]; [t C] = ode45(@KineticEqs,tspan,C0,[],k); y(:,1) = C(:,1); y(:,2:6) = C(:,2:6); f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2); f3 = y(:,3) - yexp(:,3); f4 = y(:,4) - yexp(:,4); f5 = y(:,5) - yexp(:,5); f6 = y(:,6) - yexp(:,6); f = [f1; f2; f3; f4;f5;f6]; % ------------------------------------------------------------------ function dCdt = KineticEqs(t,C,k) dCdt = ... [ ( -((k(1)+k(2)))*C(1) ) (k(1)*C(1)-k(3)*C(2)-k(4)*C(2) ) ( k(2)*C(1)-k(5)*C(3) ) ( -k(6)*C(4)+k(3)*C(2) ) ( k(4)*C(2)+k(5)*C(3)-k(7)*C(5) ) ( k(6)*C(4)+k(7)*C(5) )]; |
4楼2015-06-29 10:07:55
5楼2015-06-29 10:10:08
20