24小时热门版块排行榜    

查看: 1139  |  回复: 5
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

uklaoli

新虫 (初入文坛)

[求助] Matlab 动力学方程拟合的很差,郁闷,大侠们帮助一下。先谢了

我把程序和拟合的结果放在文件了。新手求助
回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈:小木虫仅提供交流平台,不对该内容负责。
    本内容由用户自主发布,如果其内容涉及到知识产权问题,其责任在于用户本人,如对版权有异议,请联系邮箱:xiaomuchong@tal.com
  • 附件 1 : 求助.pdf
  • 2013-06-13 20:20:42, 143.12 K

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

uklaoli

新虫 (初入文坛)

引用回帖:
5楼: Originally posted by hucanboy1 at 2013-06-14 19:28:10
楼主你好,你能不能说明一下你这个是打算仿真一个什么过程?然后你希望是什么结果?

催化裂解过程中催化剂的失活,得到K值
6楼2013-06-14 22:18:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答

uklaoli

新虫 (初入文坛)

程序是这样的
function myparafit
clear all
clc
warning off
close all
k0 = [0.0000000005  0.0000000005  0.0000000005];         % initial value
lb = [0 0 0 0];                   % ????
ub = [1 1 1 1];    % ????
x0 = [100 0  0  0 25];
yexp=[0   100   0   0      0
      20  74.1  5   3.1    17.8
      40  41.4  5   12.3   41.3
      60  7.2   7.5 20.9   64.4
      90  2.9   15  20.4   62];
%
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
k_fm= k;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
k_ls = k;

k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
k_fmls = k;


tspan=[0 20 40 60 90];

[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
plot(t,x(:,1),'*-',t,yexp(:,2),'o-');
figure;
plot(t,x(:,2:4),'*-');
hold on
plot(t,yexp(:,3:5),'o:');
figure
hold on
k_fmls;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = ObjFunc7Fmincon(k,x0,yexp)
tspan=[0 20 40 60 90];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2:5) = x(:,1:4);
f =sum((y(:,2)-yexp(:,1)).^2)+sum((y(:,3)-yexp(:,2)).^2)+sum((y(:,4)-yexp(:,3)).^2)+sum((y(:,5)-yexp(:,4)).^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = ObjFunc7LNL(k,x0,yexp)
tspan=[0 20 40 60 90];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3) = x(:,2);
y(:,4) = x(:,3);
y(:,5) = x(:,4);
f1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
f3 = y(:,4) - yexp(:,4);
f4 = y(:,5) - yexp(:,5);
f = [f1; f2; f3; f4];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dxdt = KineticEqs(t,x,k)
% global Alpha k_d
W=25;
m0=100;
k_d=0.01;
P=x(1);
G=x(2);
L=x(3);
W=x(4);
Alpha=x(5);
dxdt1=-(k(1)+k(2)+k(3))*P.*Alpha;
dxdt2=k(1)*P.*Alpha;
dxdt3=k(2)*P.*Alpha;
dxdt4=k(3)*P.*Alpha;
dxdt5=-k_d*exp(-k_d*x(3)*(100/25));
dxdt=[dxdt1;dxdt2;dxdt3;dxdt4;dxdt5];
2楼2013-06-13 23:17:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

uklaoli

新虫 (初入文坛)

请各位高人看一下了,谢谢
3楼2013-06-14 16:11:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

uklaoli

新虫 (初入文坛)

我的残差825.104
4楼2013-06-14 18:52:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-09-02 00:19 by xiangy672
[基金申请] 学科评审组评审是指会评吗? +4 瞬息宇宙 2026-08-31 4/200 2026-09-01 14:58 by jiaoxg
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +13 Kittylucky 2026-08-27 14/700 2026-09-01 11:06 by feng6531
[论文投稿] 小白求助 投论文要求的highlights应该如何写 5+3 l1963982152 2026-08-29 4/200 2026-09-01 09:04 by 北京莱茵编辑
[基金申请] 国社科又开始会评了,不知道这次命运如何 +7 雨打竹帘 2026-08-30 11/550 2026-08-31 23:16 by hittle2008
[基金申请] 基金未中,这种答复是模板吗? +6 zhaosm1982 2026-08-27 7/350 2026-08-31 21:18 by qdxxmc
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 基金不中,共勉 +12 eulota 2026-08-26 12/600 2026-08-31 08:42 by ZJTJZ
[基金申请] 为什么到现在没收到通知? +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[考博] 找导师 +6 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +15 晴天加油 2026-08-26 16/800 2026-08-29 18:28 by symmetry
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
信息提示
请填处理意见