应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 专业学科区 » 数学 » 计算数学 » 【求助】用matlab最优化方法进行参数拟合
51/1返回列表
查看: 3643  |  回复: 17
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 1 个 数学EPI ，点击这里进行查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

wodaifei

银虫 (小有名气)

[交流] 【求助】用matlab最优化方法进行参数拟合已有7人参与

各位学哥、学姐，大侠们！！
    小弟最近遇到一个难题：具体为 dyA/dt=-(k1+k2+k3+k4)yA
                                                         dyB/dt=k1yA-(k5+k6+k7)yB
                                                         dyC/dt=k2yA+k5yB-(k8+k9)yC
                                                         dyD/dt=k3yA+k6yB+k8yC
                                                         dyE/dt=k4yA+k7yB+k9yC
其中y为各组分的质量分率; k为反应速率常数
其中实验数据为：反应时间        A        B        C        D        E
                     0.833       0.0632     0．699       0.187    0.0256      0.0220
                             1.25        0.055        0.722        0.199    0.0305      0.0281
                            2.5         0.0502        0.769        0.201   0.0328        0.031
根据实验数据希望能用最优化法回归出其中的九个参数（最好能用matlab），小弟刚接触这方面的知识，遇到了困难，怎么也估计不出来，希望各位学哥、学姐，大侠们能给小弟指点一下！！跪谢！！！不胜感激!!!!!1真的把小弟难住了。。
回复此楼

» 本帖已获得的红花（最新10朵）

ffjtouchlife

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
要想拥有一切，就要让自己变得足够优秀！！！
1楼2011-03-21 08:48:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

leedobb

金虫 (正式写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
怎么算出来反应常数负号，应该也可以吧
赞一下(1人) 回复此楼
有一天，我打了个瞌睡就到了这里，但我知道我掉入了时光的循环中，虽得以永生，但只有第一个循环有意义。
4楼2011-03-21 15:48:53
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 18 个回答

leedobb

金虫 (正式写手)
★ ★ ★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
nono2009(金币+5, 数学EPI+1): 鼓励详细应助 2011-03-22 07:48:24
引用回帖:
Originally posted by wodaifei at 2011-03-21 08:48:28:
各位学哥、学姐，大侠们！！
    小弟最近遇到一个难题：具体为 dyA/dt=-(k1+k2+k3+k4)yA
                                                         dyB/dt=k1yA-(k5+k6+k7)yB
                             ...

% 哦，你的题目是这个啊，
% 这只能依靠数值解，matlab肯定没有现成的程序吧。
%
% 你这样做吧。构造一个能量函数（其实就是最小二乘法啦）
%
% F =  (yA(t2)- YYA(t2,k1,k2,k3...,k9) ) ^ 2
% ..嗯，算我直接编给你吧。

%
% dyA/dt=-(k1+k2+k3+k4)yA
%                                                          dyB/dt=k1yA-(k5+k6+k7)yB
%                                                          dyC/dt=k2yA+k5yB-(k8+k9)yC
%                                                          dyD/dt=k3yA+k6yB+k8yC
%                                                          dyE/dt=k4yA+k7yB+k9yC
%   A        B        C        D        E
%                      0.833       0.0632     0．699       0.187    0.0256      0.0220
%                              1.25        0.055        0.722        0.199    0.0305      0.0281
%                             2.5         0.0502        0.769        0.201   0.0328        0.031
clear all

kk0(1:9)=[1.2,0.9,0.8, 1.3,0.6,1.6,  1.7,1.1,1.0];
t1= 1.25  - 0.833;
t2=2.5 -0.833;
y0 = [0.0632     0.699       0.187    0.0256      0.0220]';
y1 = [ 0.055        0.722        0.199    0.0305      0.0281]';
y2 = [ 0.0502        0.769        0.201   0.0328        0.031]';
dk = 0.0001;
dt = 0.1
kk=kk0
for i=1:100000
  %%% 计算 dF/dk   
for j=1:9
    %% 计算k+dk
    k=kk;k(j)=k(j)+dk;
    A = [ -(k(1)+k(2)+k(3)+k(4)) 0 0
   k(1)  -(k(5)+k(6)+k(7))  0
     k(2)  k(5)  -(k(8)+k(9)) ];
  [V,D]=eig(A);
  DD(1,1)=D(1,1); DD(2,1) = D(2,2); DD(3,1) = D(3,3);
  DD =DD;
  VD=V^-1;
  Z0 = VD*y0(1:3);
  %%% z = z0*exp( DD(i)*t)
  Z1 = Z0 .*exp(DD*t1);
  Y1 = V *Z1;
  Int_Z1 = Z1./DD-Z0./DD; %int (Z1) = Z1/DD+C
  Y1(4) = k(3)*sum(V(1,1:3)*Int_Z1)+k(6)*sum(V(2,1:3)*Int_Z1) + k(8)*sum(V(3,1:3)*Int_Z1) + y0(4);
  Y1(5) = k(4)*sum(V(1,1:3)*Int_Z1)+k(7)*sum(V(2,1:3)*Int_Z1) + k(9)*sum(V(3,1:3)*Int_Z1) + y0(4);
  Z2 = Z0 .*exp(DD*t2);
  Y2 = V *Z2;
  Int_Z2 = Z2./DD-Z0./DD; %int (Z1) = Z1/DD+C
  Y2(4) = k(3)*sum(V(1,1:3)*Int_Z2)+k(6)*sum(V(2,1:3)*Int_Z2) + k(8)*sum(V(3,1:3)*Int_Z2) + y0(4);
  Y2(5) = k(4)*sum(V(1,1:3)*Int_Z2)+k(7)*sum(V(2,1:3)*Int_Z2) + k(9)*sum(V(3,1:3)*Int_Z2) + y0(4);
  F_up = ((Y2-y2)'*(Y2-y2) + (Y1-y1)'*(Y1-y1)) /2;
  
%% 计算k=k-dk
  k=kk;k(j)=k(j)-dk;
    A = [ -(k(1)+k(2)+k(3)+k(4)) 0 0
   k(1)  -(k(5)+k(6)+k(7))  0
     k(2)  k(5)  -(k(8)+k(9)) ];
  [V,D]=eig(A);
  DD(1,1)=D(1,1); DD(2,1) = D(2,2); DD(3,1) = D(3,3);
  DD =DD;
  VD=V^-1;
  Z0 = VD*y0(1:3);
  %%% z = z0*exp( DD(i)*t)
  Z1 = Z0 .*exp(DD*t1);
  Y1 = V *Z1;
  Int_Z1 = Z1./DD-Z0./DD; %int (Z1) = Z1/DD+C
  Y1(4) = k(3)*sum(V(1,1:3)*Int_Z1)+k(6)*sum(V(2,1:3)*Int_Z1) + k(8)*sum(V(3,1:3)*Int_Z1) + y0(4);
  Y1(5) = k(4)*sum(V(1,1:3)*Int_Z1)+k(7)*sum(V(2,1:3)*Int_Z1) + k(9)*sum(V(3,1:3)*Int_Z1) + y0(4);
  Z2 = Z0 .*exp(DD*t2);
  Y2 = V *Z2;
  Int_Z2 = Z2./DD-Z0./DD; %int (Z1) = Z1/DD+C
  Y2(4) = k(3)*sum(V(1,1:3)*Int_Z2)+k(6)*sum(V(2,1:3)*Int_Z2) + k(8)*sum(V(3,1:3)*Int_Z2) + y0(4);
  Y2(5) = k(4)*sum(V(1,1:3)*Int_Z2)+k(7)*sum(V(2,1:3)*Int_Z2) + k(9)*sum(V(3,1:3)*Int_Z2) + y0(4);
  F_down = ((Y2-y2)'*(Y2-y2) + (Y1-y1)'*(Y1-y1)) /2;
  FF= (F_down+F_up)/2;
  DF(j)= (F_up-F_down)/dk/2;
end

for j=1:9
   kk(j)= kk(j) - DF(j) *dt;
end

if(mod(i,1000)==1)
   'deviation'
   FF/10
   'k1->k9', kk
end
  
  
  
  
  
  
  
end

[ Last edited by leedobb on 2011-3-21 at 15:46 ]
赞一下(6人) 回复此楼
有一天，我打了个瞌睡就到了这里，但我知道我掉入了时光的循环中，虽得以永生，但只有第一个循环有意义。
2楼2011-03-21 15:36:02
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

leedobb

金虫 (正式写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
引用回帖:
Originally posted by wodaifei at 2011-03-21 08:48:28:
各位学哥、学姐，大侠们！！
    小弟最近遇到一个难题：具体为 dyA/dt=-(k1+k2+k3+k4)yA
                                                         dyB/dt=k1yA-(k5+k6+k7)yB
                             ...

这是算出来比较好的结果
k = 0.2974    0.3073   -0.2075    0.0272    0.1690   -0.1148   -0.0798    0.4677    0.2665
此时
Y(t2)=  0.0529    0.7138    0.1869    0.0232    0.0235
y(t3)= 0.0311    0.7527    0.1864    0.0163    0.0138
t1用来当初始条件了。

唉，看起来，我太无聊了。

[ Last edited by leedobb on 2011-3-21 at 15:45 ]
赞一下(1人) 回复此楼
有一天，我打了个瞌睡就到了这里，但我知道我掉入了时光的循环中，虽得以永生，但只有第一个循环有意义。
3楼2011-03-21 15:44:50
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wodaifei

银虫 (小有名气)

哥，你太给力了！！！！！！！！

哥，顶死你！！！你太给力了。。。。。虽然动力学常数存在负值这是不对的，但问题不是你的程序，是我的实验数据。小弟跪谢！！！！！！！！！！！
赞一下 回复此楼
要想拥有一切，就要让自己变得足够优秀！！！
5楼2011-03-21 17:48:04
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 18 个回答
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭