版塊導(dǎo)航: 正在加載中...

登錄注冊

應(yīng)《網(wǎng)絡(luò)安全法》要求，自2017年10月1日起，未進(jìn)行實(shí)名認(rèn)證將不得使用互聯(lián)網(wǎng)跟帖服務(wù)。為保障您的帳號能夠正常使用，請盡快對帳號進(jìn)行手機(jī)號驗(yàn)證，感謝您的理解與支持！

24小時(shí)熱門版塊排行榜

返回列表

kidingkids

新蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 2.5
散金: 10
帖子: 5
在線: 4.1小時(shí)
蟲號: 1292152
注冊: 2011-05-11
性別: GG
專業(yè): 化學(xué)反應(yīng)工程

[求助] matlab最小二乘插值擬合，跪求大神��！

function benan5
clear all;
load cc  %load experimental data

[y_row,y_col]=size(y);
[x_row,x_col]=size(x);

global b

beta0=[2 60 1 45 1.5 90 5 45 1 100 1 0 1 1 1 1 1 1 0]
y0=[0 0 0];
lb=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
ub=[+inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf];

[beta,resnorm,residual,exitflag,output,lambda,jacobian]=...
lsqnonlin(@func,beta0,lb,ub,[],x,y0,y)
ci=nlparci(beta,residual,jacobian);

namex='W/Fa0 , g.h/mol';

for i=1:1:3
T1=180+273;
T2=200+273;
T3=220+273;
T4=240+273;
x1=x(1:8,i);
x2=x(9:16,i);
x3=x(17:24,i);
x4=x(25:32,i);
b=3*i-2;
x5=linspace(0,500,1000);
tspan1=[0 max(x1)];
tspan2=[0 max(x2)];
tspan3=[0 max(x3)];
tspan4=[0 max(x4)];

y1=y(1:8,b);
y2=y(9:16,b);
y3=y(17:24,b);
y4=y(25:32,b);
y5=y(1:8,b+1);
y6=y(9:16,b+1);
y7=y(17:24,b+1);
y8=y(25:32,b+1);
y9=y(1:8,b+2);
y10=y(9:16,b+2);
y11=y(17:24,b+2);
y12=y(25:32,b+2);

[t1 yy1]=ode23(@modeleqs,tspan1,y0,[],beta,T1);
[t2 yy2]=ode23(@modeleqs,tspan2,y0,[],beta,T2);
[t3 yy3]=ode23(@modeleqs,tspan3,y0,[],beta,T3);
[t4 yy4]=ode23(@modeleqs,tspan4,y0,[],beta,T4);

yc1=spline(t1,yy1(:,1),x1);
yc2=spline(t2,yy2(:,1),x2);
yc3=spline(t3,yy3(:,1),x3);
yc4=spline(t4,yy4(:,1),x4);

yc5=spline(t1,yy1(:,2),x1);
yc6=spline(t2,yy2(:,2),x2);
yc7=spline(t3,yy3(:,2),x3);
yc8=spline(t4,yy4(:,2),x4);

yc9=spline(t1,yy1(:,3),x5);
yc10=spline(t2,yy2(:,3),x5);
yc11=spline(t3,yy3(:,3),x5);
yc12=spline(t4,yy4(:,3),x5);

figure(6*i-5)
plot(x1,yc1,'r-',x1,y1,'o');hold on;
plot(x2,yc2,'b-',x2,y2,'s');hold on;
plot(x3,yc3,'m-',x3,y3,'v');hold on;
plot(x4,yc4,'g-',x4,y4,'x');hold off;
xlabel(namex)
ylabel('苯胺轉(zhuǎn)化率')

figure(6*i-4)
plot(x1,yc5,'r-',x1,y5,'o');hold on;
plot(x2,yc6,'b-',x2,y6,'s');hold on;
plot(x3,yc7,'m-',x3,y7,'v');hold on;
plot(x4,yc8,'g-',x4,y8,'x');hold off;
xlabel(namex)
ylabel('二環(huán)己胺收率')

figure(6*i-3)
plot(x5,yc9,'r-',x1,y9,'o');hold on;
plot(x5,yc10,'b-',x2,y10,'s');hold on;
plot(x5,yc11,'m-',x3,y11,'v');hold on;
plot(x5,yc12,'g-',x4,y12,'x');hold off;
xlabel(namex)
ylabel('N-環(huán)己基苯胺收率')
%{
figure(6*i-2)
b=3*i-2
d1=residual(96*(i-1)+1:96*(i-1)+8,1);hold on;
d2=residual(96*(i-1)+9:96*(i-1)+16,1);hold on;
d3=residual(96*(i-1)+17:96*(i-1)+24,1);hold on;
d4=residual(96*(i-1)+25:96*(i-1)+32,1);hold off;
plot(y1,d1,'ko',y2,d2,'kv',y3,d3,'k*',y4,d4,'ks');
namex1='苯胺轉(zhuǎn)化率';
namey1(1,1:length('residual values'))='Residual values';
xlabel(namex1)
ylabel(namey1)

figure(6*i-1)
d5=residual(96*(i-1)+33:96*(i-1)+40,1);hold on;
d6=residual(96*(i-1)+32+9:96*(i-1)+32+16,1);hold on;
d7=residual(96*(i-1)+32+17:96*(i-1)+32+24,1);hold on;
d8=residual(96*(i-1)+32+25:96*(i-1)+32+32,1);hold off;
plot(y5,d5,'ko',y6,d6,'kv',y7,d7,'k*',y8,d8,'ks');
namex2='二環(huán)己胺收率';
namey2(1,1:length('residual values'))='Residual values';
xlabel(namex2)
ylabel(namey2)

figure(6*i)
d9=residual(96*(i-1)+64+1:96*(i-1)+64+8,1);hold on;
d10=residual(96*(i-1)+64+9:96*(i-1)+64+16,1);hold on;
d11=residual(96*(i-1)+64+17:96*(i-1)+64+24,1);hold on;
d12=residual(96*(i-1)+64+25:96*(i-1)+64+32,1);hold off;
plot(y9,d9,'ko',y10,d10,'kv',y11,d11,'k*',y12,d12,'ks');
namex3='N-環(huán)己基苯胺收率';
namey3(1,1:length('residual values'))='Residual values';
xlabel(namex3)
ylabel(namey3)
%}
end

%output of parameters needed
fprintf('\n estimated parameters by lsqnonlin():\n')
fprintf('\t k1=%10.6f±%10.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\t k2=%10.6f±%10.4f\n',beta(2),ci(2,2)-beta(2))
fprintf('\t k3=%10.6f±%10.4f\n',beta(3),ci(3,2)-beta(3))
fprintf('\t k4=%10.6f±%10.4f\n',beta(4),ci(4,2)-beta(4))
fprintf('\t k5=%10.6f±%10.4f\n',beta(5),ci(5,2)-beta(5))
fprintf('\t k6=%10.4f±%10.4f\n',beta(6),ci(6,2)-beta(6))
fprintf('\t k7=%10.4f±%10.4f\n',beta(7),ci(7,2)-beta(7))
fprintf('\t k8=%10.4f±%10.4f\n',beta(8),ci(8,2)-beta(8))
fprintf('\t k9=%10.4f±%10.4f\n',beta(9),ci(9,2)-beta(9))
fprintf('\t k10=%10.4f±%10.4f\n',beta(10),ci(10,2)-beta(10))
fprintf('\t k11=%10.6f±%10.4f\n',beta(11),ci(11,2)-beta(11))
fprintf('\t k12=%10.6f±%10.4f\n',beta(12),ci(12,2)-beta(12))
fprintf('\t k13=%10.6f±%10.4f\n',beta(13),ci(13,2)-beta(13))
fprintf('\t k14=%10.6f±%10.4f\n',beta(14),ci(14,2)-beta(14))
fprintf('\t k15=%10.6f±%10.4f\n',beta(15),ci(15,2)-beta(15))
fprintf('\t k16=%10.4f±%10.4f\n',beta(16),ci(16,2)-beta(16))
fprintf('\t k17=%10.4f±%10.4f\n',beta(17),ci(17,2)-beta(17))
fprintf('\t k18=%10.4f±%10.4f\n',beta(18),ci(18,2)-beta(18))
fprintf('\t k19=%10.4f±%10.4f\n',beta(19),ci(19,2)-beta(19))
fprintf(' the sum of the residual squares is:%.4e\n\n',sum((residual).^2))

%=======================================================
function f=func(beta,x,y0,y)
%define objective function
global b
for i=1:1:3
T1=180+273;
T2=200+273;
T3=220+273;
T4=240+273;
x1=x(1:8,i);
x2=x(9:16,i);
x3=x(17:24,i);
x4=x(25:32,i);
b=3*i-2;

tspan1=[0 max(x1)];
tspan2=[0 max(x2)];
tspan3=[0 max(x3)];
tspan4=[0 max(x4)];

[t1 yy1]=ode23(@modeleqs,tspan1,y0,[],beta,T1);
[t2 yy2]=ode23(@modeleqs,tspan2,y0,[],beta,T2);
[t3 yy3]=ode23(@modeleqs,tspan3,y0,[],beta,T3);
[t4 yy4]=ode23(@modeleqs,tspan4,y0,[],beta,T4);

y1=y(1:8,b);
y2=y(9:16,b);
y3=y(17:24,b);
y4=y(25:32,b);

y5=y(1:8,b+1);
y6=y(9:16,b+1);
y7=y(17:24,b+1);
y8=y(25:32,b+1);

y9=y(1:8,b+2);
y10=y(9:16,b+2);
y11=y(17:24,b+2);
y12=y(25:32,b+2);

j=1;
yc1=spline(t1,yy1(:,j),x1);
yc2=spline(t2,yy2(:,j),x2);
yc3=spline(t3,yy3(:,j),x3);
yc4=spline(t4,yy4(:,j),x4);
%苯胺轉(zhuǎn)化率擬合
j=2;
yc5=spline(t1,yy1(:,j),x1);
yc6=spline(t2,yy2(:,j),x2);
yc7=spline(t3,yy3(:,j),x3);
yc8=spline(t4,yy4(:,j),x4);
   %二環(huán)己胺收率擬合
j=3;
yc9=spline(t1,yy1(:,j),x1);
yc10=spline(t2,yy2(:,j),x2);
yc11=spline(t3,yy3(:,j),x3);
yc12=spline(t4,yy4(:,j),x4);
%N-環(huán)己基苯胺收率擬合

f11=y1-yc1;
f12=y2-yc2;
f13=y3-yc3;
f14=y4-yc4;
f21=y5-yc5;
f22=y6-yc6;
f23=y7-yc7;
f24=y8-yc8;
f31=y9-yc9;
f32=y10-yc10;
f33=y11-yc11;
f34=y12-yc12;

g(i,

=[f11' f12' f13' f14' f21' f22' f23' f24' f31' f32' f33' f34'];
end
f=[g(1,

';g(2,

';g(3,

'];
%===================================================================================
function dydt=modeleqs(t,yt,beta,T)
global b
P=0.3;
N=5*(b+2)/3+5;%N為氫胺比
X1=yt(1);
Y4=yt(2);
Y6=yt(3);
X1
Y4
Y6
P1=P.*(1-X1)./(N+1-3.*X1+3.*Y6);
P2=P.*(N-3.*X1+3.*Y6)./(N+1-3.*X1+3.*Y6);
P3=P.*(X1-Y4-2.*Y6)./(N+1-3.*X1+3.*Y6);
P4=P.*Y4./(2.*N+2-6.*X1+6.*Y6);
P5=P.*(Y4+2.*Y6)./(2.*N+2-6.*X1+6.*Y6);
P6=P.*Y6./(N+1-3.*X1+3.*Y6);
P1
P3
P4
P5
P6
K(1)=beta(1).*10.^6.*exp(-beta(2).*10.^3./(8.314.*T)); %k(1)、k(2)第一步反應(yīng)的指前因子、表觀活化能
K(2)=beta(3).*10.^7.*exp(-beta(4).*10.^3./(8.314.*T)); %k(3)、k(4)第二步反應(yīng)的指前因子、表觀活化能
K(3)=beta(5).*10.^10.*exp(-beta(6).*10.^3./(8.314.*T)); %k(5)、k(6)第三步反應(yīng)的指前因子、表觀活化能
K(4)=beta(7).*10.^11.*exp(-beta(8).*10.^3./(8.314.*T)); %k(7)、k(8)第四步反應(yīng)的指前因子、表觀活化能
K(5)=beta(9).*10.^12.*exp(-beta(10).*10.^3./(8.314.*T)); %k(9)、k(10)第五步反應(yīng)的指前因子、表觀活化能
r1=K(1).*P1.^beta(11).*P2.^beta(12);
r2=K(2).*P3.^beta(13)-K(3).*P4.^beta(14).*P5.^beta(15);
r3=K(4).*P1.^beta(16).*P3.^beta(17);
r4=K(5).*P6.^beta(18).*P2.^beta(19);
R1=r1+r3
R4=r2+r4
R6=r3-r4
dydt=[r1+r3;r2+r4;r3-r4];
求大神檢查一下程序是否有問題。我現(xiàn)在擬合不出來結(jié)果，急！實(shí)驗(yàn)數(shù)據(jù)在附件里

回復(fù)此樓

» 本帖附件資源列表

歡迎監(jiān)督和反饋：小木蟲僅提供交流平臺，不對該內(nèi)容負(fù)責(zé)。
本內(nèi)容由用戶自主發(fā)布，如果其內(nèi)容涉及到知識產(chǎn)權(quán)問題，其責(zé)任在于用戶本人，如對版權(quán)有異議，請聯(lián)系郵箱：xiaomuchong@tal.com
附件 1 : cc.mat

2013-11-06 10:16:48, 2.17 K

» 猜你喜歡

光學(xué)工程學(xué)碩調(diào)劑信息已經(jīng)有32人回復(fù)
歡迎加入課題組已經(jīng)有0人回復(fù)
物理學(xué)I論文潤色/翻譯怎么收費(fèi)? 已經(jīng)有233人回復(fù)
散金幣，求好運(yùn)，祝面上順利！已經(jīng)有79人回復(fù)
華東師范大學(xué)芯片設(shè)計(jì)徐瓏真實(shí)水平咋樣已經(jīng)有1人回復(fù)
調(diào)劑已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項(xiàng)目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有0人回復(fù)
中德博士后2026項(xiàng)目5月30日截止- 國內(nèi)博士后入站、德國兩年全職已經(jīng)有0人回復(fù)
面向全球申請人的中德博士后2026項(xiàng)目5月30日截止- 國內(nèi)博后入站德國兩年全職已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項(xiàng)目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有2人回復(fù)
重慶交大材料學(xué)院26年碩士招生復(fù)試實(shí)施細(xì)則已出。歡迎加入光子學(xué)微結(jié)構(gòu)與器件課題組。已經(jīng)有0人回復(fù)

1樓 2013-11-06 10:22:11

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

智能機(jī)器人

Robot (super robot)

我們都愛小木蟲

找到一些相關(guān)的精華帖子，希望有用哦~

matlab插值擬合已經(jīng)有8人回復(fù)
概率論和數(shù)理統(tǒng)計(jì)的Matlab實(shí)現(xiàn) 已經(jīng)有261人回復(fù)
求助Matlab最小二乘法非線性擬合已經(jīng)有4人回復(fù)
部分函數(shù)曲線和幾個(gè)點(diǎn)怎么擬合已經(jīng)有9人回復(fù)
基于MATLAB Simulink的系統(tǒng)仿真技術(shù)與應(yīng)用--薛定宇--PDF 已經(jīng)有149人回復(fù)
求助，一自變量二因變量擬合matlab該怎么實(shí)現(xiàn)呢？已經(jīng)有9人回復(fù)
matlab對數(shù)據(jù)進(jìn)行最小二乘擬合，若果數(shù)據(jù)為負(fù)數(shù)怎么擬合？急急急！已經(jīng)有3人回復(fù)
MATLAB微分方程參數(shù)擬合問題，求大神已經(jīng)有7人回復(fù)
GS+和Arcgis變異函數(shù)擬合已經(jīng)有9人回復(fù)
基于最小二乘互相關(guān)算法的圖像定位匹配研究（論文討論）已經(jīng)有8人回復(fù)
非線性數(shù)據(jù)擬合的數(shù)學(xué)原理！�。。。�？？？？已經(jīng)有30人回復(fù)
matlab如何編寫共享參數(shù)擬合程序已經(jīng)有8人回復(fù)
擬合和回歸有什么區(qū)別已經(jīng)有8人回復(fù)
樣條插值與樣條擬合的區(qū)別已經(jīng)有5人回復(fù)
50金幣求一MATLAB的擬合程序已經(jīng)有11人回復(fù)
請問將一批離散的點(diǎn)(二維平面上)擬合為曲線目前有哪些比較新的方法? 已經(jīng)有3人回復(fù)
【求助】Matlab非線性最小二乘擬合活度系數(shù)模型（Willson、NRTL、UNIQUAC）已經(jīng)有17人回復(fù)

點(diǎn)擊這里搜索更多相關(guān)資源

科研從小木蟲開始，人人為我，我為人人

相關(guān)版塊跳轉(zhuǎn) 我要訂閱樓主 kidingkids 的主題更新