24小时热门版块排行榜    

查看: 583  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 时间戳变了,能看出什么问题? +18 基诺咪客 2026-08-17 23/1150 2026-08-20 17:19 by Godzela
[基金申请] 只有每年这种时候来逛逛小木虫 +21 yaoyewhu2008 2026-08-20 22/1100 2026-08-20 17:03 by yaoyewhu2008
[基金申请] 今天系统多次维护,明天很可能放榜! +9 zju2000 2026-08-16 10/500 2026-08-20 15:56 by Leogzhya
[基金申请] 放榜前的不淡定 40+3 snowwithsea 2026-08-19 12/600 2026-08-20 15:56 by lfy8008
[基金申请] 言之凿凿,明天官网会通知放榜? +14 医学老男孩 2026-08-20 16/800 2026-08-20 15:47 by LNP@mRNA
[基金申请] 今天基金会出结果吗?20260819 +14 kkkl_v 2026-08-19 15/750 2026-08-20 08:21 by sunzitan
[基金申请] filecode,4个jtjc了 +11 ziyangfang 2026-08-19 13/650 2026-08-19 21:28 by aasahr
[教师之家] 为什么余额宝的年化利率越来越低?主要原因有哪些? +5 瞬息宇宙 2026-08-15 5/250 2026-08-19 20:23 by super2002521
[基金申请] 今天放榜没戏了吧 +4 yuleib84 2026-08-19 5/250 2026-08-19 19:52 by hhs666
[基金申请] 今天放榜吗? +14 布布和一二 2026-08-19 15/750 2026-08-19 18:07 by gltch
[基金申请] filecode=后面第一个是大写字母 +10 wangze12014 2026-08-14 12/600 2026-08-19 16:56 by 苦难博士
[基金申请] 2027广东省杰青 +4 奶牛小黑 2026-08-15 10/500 2026-08-19 11:02 by wanfengnew
[基金申请] 朋友圈看到的 +6 wangzilk 2026-08-18 8/400 2026-08-19 10:55 by Haru815
[基金申请] 快农历七夕节了,轻松一下,男人悄悄话,女施主请不要进来。 +6 Tide man 2026-08-14 7/350 2026-08-19 09:56 by ZJTJZ
[基金申请] 什么时候开奖? +6 CrisMessi 2026-08-18 6/300 2026-08-19 08:06 by Equinoxhua
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
[基金申请] 感觉是下周放榜了 +6 angus9576 2026-08-17 11/550 2026-08-17 23:57 by angus9576
[精细化工] 招聘 金属平磨液,抛光液研发工程师 +3 小天0311 2026-08-14 3/150 2026-08-16 07:31 by H9PLUS
信息提示
请填处理意见