| 查看: 960 | 回复: 1 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
_romantic_镇木虫 (著名写手)
|
[求助]
初学者求助，matlab模仿程序已有1人参与
|
|
|
function KineticsEst5 % 动力学ODE方程模型的参数估计 clear all clc k0 = [0 0 0 0]; % 参数初值 lb = [ -inf -inf -inf -inf ]; % 参数下限 ub = [+inf +inf +inf +inf ]; % 参数上限 x0 = [8.5 28.8 27.6]; tspan=[0 0.222 0.333 0.444]; yexp = [8.5000 28.8000 27.6000 4.8000 23.2000 35.3000 4.2000 21.6000 36.5000 4.0000 21.2000 37.3000]; % yexp: 实验数据[x1 x2 x3] % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,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(' 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,[],x0,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,[],x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n') Output % ------------------------------------------------------------------ function f = ObjFunc4Fmincon(k,x0,yexp) tspan = [0.00 : 0.01 : 0.44]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1:3) = x(:,1:3); f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2) ... + sum((y(:,3)-yexp(:,3)).^2) ; % 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0 = k_fmincon; [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n') Output % ------------------------------------------------------------------ function f = ObjFunc4LNL(k,x0,yexp) tspan = [0.00 : 0.01 : 0.44]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1:3) = x(:,1:3); f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2); f3 = y(:,3) - yexp(:,3); f = [f1; f2; f3]; % ------------------------------------------------------------------ function dxdt = KineticEqs(t,x,k) dxdt = ... [ (k(1)*x(2)- k(1)*x(1)) (k(1)*x(2)+(k(3)-k(2)-k(4))*x(3)) (k(3)*x(2)-k(4)*x(3)) ]; |
回复此楼» 猜你喜欢
求助将cif文件转化为jade
已经有1人回复
-
金属氧化物薄膜空气退火和真空退火对表面粗糙度的影响
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有119人回复
-
Apex2单晶测试
已经有0人回复
-
结构建模
已经有0人回复
-
声子谱计算问题
已经有2人回复
-
已结束
已经有23人回复
-
Error in routine lambda (100): wrong # or too many modes
已经有0人回复
-
祈福国自然面上
已经有55人回复
-
请问castep计算出的能带图横轴数据缩放因子怎么计算
已经有1人回复
200