24小时热门版块排行榜    

查看: 531  |  回复: 0

2011207156

金虫 (小有名气)

[求助] matlab拟合动力学参数,非常谢谢~

function butongshike201506171
clear all;clc
format long
tspan=(0:10:390);
yexp=[0.046125761496241
   0.078139373146412
   0.117158327134899
   0.159846402230616
   0.199412869433285
   0.231111460229405
   0.259814231746068
   0.282697037233637
   0.306493661349134
   0.326396561089024
   0.345295684999190
   0.360444423309371
   0.377418750129898
   0.395250807577438
   0.414965633800510
   0.430017989495138
   0.445923212579888
   0.457995753035510
   0.469186259425871
   0.480284344291621
   0.491524402631167
   0.501658351934599
   0.510935126510884
   0.520202070112304
   0.531309455870521
   0.537782255717000
   0.547944130180314
   0.558094213007469
   0.565468701311900
   0.573757577102420
   0.581007002009939
   0.589392949772071
   0.597541754780127
   0.607751809867969
   0.613129267171725
   0.620580135887054
   0.626072239599544
   0.632475341850705
   0.638742438783615];
beta0=[67971 0.000839];
y0=0;
lb=[-1 -1]*1e6;
ub=[1e9 1e9];   
yy=[0 yexp'];
k0=beta0;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\t待拟合参数 k1 = %.6f\n',k(1))
fprintf('\t待拟合参数 k2= %.6f\n',k(2))
fprintf('  残差平方和= %.6f\n\n',resnorm)
ts=0:10:max(tspan);
[ts ys]=ode45(@KineticsEqs,ts,y0,[],k)
[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k);
y=XXsim(2:end);
xexp=yexp;
R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);
fprintf('\n\t决定系数R-Square = %.6f',R2);
figure(1)
plot(ts,ys,'b',tspan,yy,'or'),legend('计算值','实验值','Location','best');
xlabel('x'),ylabel('y')


%---------------------------------------------------------
function f = ObjFunc(k,tspan,y0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;
ysim = Xsim(2:end);
size(ysim);
size(yexp);
f=ysim-yexp;
%----------------------------------------------------------

function dydt = KineticsEqs(t,y,k)
beta(1)=k(1);
beta(2)=k(2);
tspan=(0:10:390);
h=0.068031;
p=[3.472000000000000
   3.519000000000000
   3.548000000000000
   3.589000000000000
   3.634000000000000
   3.676000000000000
   3.710000000000000
   3.740000000000000
   3.765000000000000
   3.789000000000000
   3.811000000000000
   3.832000000000000
   3.847000000000000
   3.865000000000000
   3.885000000000000
   3.906000000000000
   3.921000000000000
   3.938000000000000
   3.952000000000000
   3.964000000000000
   3.977000000000000
   3.988000000000000
   4.000000000000000
   4.010000000000000
   4.020000000000000
   4.032000000000000
   4.039000000000000
   4.050000000000000
   4.061000000000000
   4.069000000000000
   4.078000000000000
   4.087000000000000
   4.095000000000000
   4.105000000000000
   4.115000000000000
   4.122000000000000
   4.129000000000000
   4.135000000000000
   4.142000000000000
   4.150000000000000]';
D=p*h;
sp=csaps(tspan,D,1);
pp=fnder(sp);
dD=fnval(pp,tspan);
dydt= (83.38-0.0014*exp(1600*(1/277.55-1/298.15))...
*interp1(tspan,p,t,'pchip')*10000-1/beta(2)...
*interp1(tspan,D,t,'pchip')*4.229*(0.0598...
-(y*0.0150*(54.2+17.68)+(0.236204-interp1(tspan,D,t,'pchip'))...
*0.0598*4.229)/(45.07-interp1(tspan,D,t,'pchip')*4.229))...
*interp1(tspan,dD,t,'pchip')/(-45.07+interp1(tspan,D,t,'pchip')...
*4.229)/0.0150-(0.0598-(y*0.0150*(54.2+17.68)+(0.236204-...
interp1(tspan,D,t,'pchip'))*0.0598*4.229)/...
(45.07-interp1(tspan,D,t,'pchip')*4.229))*...
interp1(tspan,dD,t,'pchip')/0.0150/beta(1))/...
(1/beta(1)+1/beta(2)-1/beta(2)*interp1(tspan,D,t,'pchip')...
*71.88/(-45.07+4.229*interp1(tspan,D,t,'pchip')));
以上是求动力学参数beta(1)和beta(2)的代码,p,D,和dD是随时间变化的已知量,不知以上程序运行时interp1(tspan,D,t,'pchip')能否对应上不同时刻的D。非常谢谢。
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子,希望有用哦~

科研从小木虫开始,人人为我,我为人人
相关版块跳转 我要订阅楼主 2011207156 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 求调剂,一志愿:南京航空航天大学大学 ,080500材料科学与工程学硕,总分289分 +4 @taotao 2026-03-19 4/200 2026-03-20 22:14 by JourneyLucky
[考研] 085600材料与化工 +8 安全上岸! 2026-03-16 8/400 2026-03-20 22:13 by luoyongfeng
[考研] 南昌大学材料专硕311分求调剂 +3 77chaselx 2026-03-20 3/150 2026-03-20 22:03 by luoyongfeng
[考研] 281求调剂(0805) +14 烟汐忆海 2026-03-16 25/1250 2026-03-20 15:47 by yuncha
[考研] 0856调剂,是学校就去 +8 sllhht 2026-03-19 9/450 2026-03-20 14:25 by 无懈可击111
[考研] 求调剂 +3 暗涌afhb 2026-03-16 3/150 2026-03-20 00:28 by 河南大学校友
[考研] 081700化工学硕调剂 +3 【1】 2026-03-16 3/150 2026-03-19 23:40 by edmund7
[考研] 320求调剂0856 +3 不想起名字112 2026-03-19 3/150 2026-03-19 22:53 by 学员8dgXkO
[考研] 085600材料与化工求调剂 +6 绪幸与子 2026-03-17 6/300 2026-03-19 13:27 by houyaoxu
[考研] 0817调剂 +3 没有答案_ 2026-03-14 3/150 2026-03-19 09:51 by Xu de nuo
[考研] 297求调剂 +8 戏精丹丹丹 2026-03-17 8/400 2026-03-18 14:30 by laoshidan
[考研] 301求调剂 +9 yy要上岸呀 2026-03-17 9/450 2026-03-18 08:58 by 无际的草原
[硕博家园] 湖北工业大学 生命科学与健康学院-课题组招收2026级食品/生物方向硕士 +3 1喜春8 2026-03-17 5/250 2026-03-17 17:18 by ber川cool子
[考研] 材料专硕326求调剂 +6 墨煜姒莘 2026-03-15 7/350 2026-03-17 17:10 by ruiyingmiao
[考研] 一志愿,福州大学材料专硕339分求调剂 +3 木子momo青争 2026-03-15 3/150 2026-03-17 07:52 by laoshidan
[考研] 11408 一志愿西电,277分求调剂 +3 zhouzhen654 2026-03-16 3/150 2026-03-17 07:03 by laoshidan
[考研] 304求调剂 +5 素年祭语 2026-03-15 5/250 2026-03-16 17:00 by 我的船我的海
[考研] 326求调剂 +3 mlpqaz03 2026-03-15 3/150 2026-03-16 07:33 by Iveryant
[考研] 085601材料工程315分求调剂 +3 yang_0104 2026-03-15 3/150 2026-03-15 10:58 by peike
[考研] 080500,材料学硕302分求调剂学校 +4 初识可乐 2026-03-14 5/250 2026-03-14 21:08 by peike
信息提示
请填处理意见