?上一篇文章：
?下一篇文章：
& &评论摘要(共 0 条，得分 0 分，平均 0 分)小木虫 --- 500万硕博科研人员喜爱的学术科研平台
&&查看话题
Matlab 对 ODE的参数进行 拟合 求助@月只蓝
运行 一直提示 Not enough input arguments 以下是我的代码
1 Call lsqnonlin
tic % time check
%Experimental data
tspan = ';
& &% Concentration of Mg2+&&mol/L
% initial parameter identification
kn = 10^2.6; %Crystallization kinetics
kg = 8.06e-2;
X_1=1; % mass transfer coefficient enhancement factor
teta_0 = ; %just collect in 1 vector sampling instants t and measured B
%initial values for ODE
% Call the optimizer:
& &= lsqnonlin(@Objfun_ffMgOH2_CO2_H2O_with_precipitation,teta_0,lb,ub,,tspan,x0,yexp);
ci = nlparci(teta_opt,residual,jacobian);
kn = teta_opt(1);
kg = teta_opt(2);
n&&= teta_opt(3);
g&&= teta_opt(4);
X_1 = teta_opt(5);
X_2 = teta_opt(6);
X_3 = teta_opt(7);
%Optimized solution
= ode15s(@ode_ffMgOH2_CO2_H2O_with_precipitation,tsim,x0,,kn,kg,n,g,X_1,X_2,X_3);
plot(tspan,y_1,'o',tsim,x(:,7)+x(:,9)) %plot the data vs solution
title('Parameter fitting of Mg2+ in the bulk,mol/L');
PH=-log10(x(:,4));
plot(tspan,y_2,'o',tsim,PH) %plot the data vs solution
title('Parameter fitting of PH');
toc %time check
2 Call Objfun
%Call 'Objfun' function:
function lsq = Objfun_ffMgOH2_CO2_H2O_with_precipitation(teta,tspan,x0,yexp)
= ode15s(@ode_ffMgOH2_CO2_H2O_with_precipitation,tspan,x0,,teta);
%Concentration of Mg2+
y_cal(:,1) = x(:,7)+x(:,9);
y_cal(:,2) = -log10(x(:,4));
3 Call ode_function
% the reaction of Calc hydroxide and Carbon dioxide
function dx=ode_ffMgOH2_CO2_H2O_with_precipitation(t,x,kn,kg,n,g,X_1,X_2,X_3)
dx = zeros(13,1);
% x(1) is the concentration of liquid CO2(l)
% initial value = 0 mol/L
% x(2) is the concentration of gas CO2(g)
% initial value = 0 mol/L
% x(3) is the concentration of Hydroxy&&in bulk
% initial value = 3.3e-4 mol/L
% x(4) is the concentration of Hydrogen&&in bulk
% initial value = 3.03e-11 mol/L
% x(5) is the concentration of Bicarbonate&&
% initial value = 0 mol/L
% x(6) is the concentration of carbonate&&
% initial value = 0 mol/L
% x(7) is the concentration of Magnesium
% initial value = 1.65e-4&&mol/L
% x(8) is the mass of Mg(OH)2
% initial value = 100&&g
% x(9) is the concentraiton of liquid MgCO3
% initial value = 0&&mol/L
% x(10) is o monent which is number density of crystal n(L)
% initial value = 0
% x(11) is 1 monent which is number density of crystal n(L)
% initial value = 0
% x(12) is 2 monent which is number density of crystal n(L)
% initial value = 0
% x(13) is 3 monent which is number density of crystal n(L)
% initial value = 0
Qg = 1.67e-2;& && &&&% L/s& && && & 1L/min
N&&= 650;&&
vl = 2.826;& && && & % vl is the dispersion volum& && && && && && &L
vg = 2.047*10^(-5)*N^(0.579)*Qg^0.63*&&
& && && && && && && &% vg is the gas volum in the gas-liquid mixture m3
& && && &&&%agitating speed rpm
P = 1;& && && && && &% P is the atmospheric pressure& && && && && & atm
pi = 3.14;
kla = 0.0304;& && &&&% kla is the mass transfer coefficient of CO2&&s^(-1)
& && && && && && && &%(9.6218e-5)*N^0.4316*Qg^0.3840*c^0.1117*E;&&
kH = 0.035;& && && & % kh is the Henry's constant& && && && && &mol/(L*atm)
M_CO2 = 44;& && && & % M_CO2 is the co2 molar mass& && && && &&&g/mol
M_MgOH = 58;& && && &% molar mass of Mg(OH)2& && && && && &g/mol
M_MgCO3 = 84;& && &&&% molar mass of M_MgCO3& & g/mol
dd_CO2 = 1.8;& && &&&% dd_CO2 is the density of CO2
dd_MgOH = 2.34e3;& & % dd_Mg(OH)2 is the density of Mg(OH)2& && && & g/L
dd_MgCO3 = 2.96e3;& &% dd_Mg(OH)2 is the density of MgCO3& && && & g/L
dd_mol_MgCO3 = dd_MgCO3/M_MgCO3; % mol/L
ddslurry = 1e3;& && &% density of slurry g/L& && && && && && && &% ddCO2 = 1960; ddCO2 is the density of CO2
c = 1;& && && && && &% volume fraction of CO2 in the inlet gas& && && && && && && &
CO2in = dd_CO2/M_CO2*c;% CO2in is the concentration of CO2 entering the reactor mol/L
& && && && && && && &% CO2in = ddCO2/MCO2*0.15 =.05 =6.75&&mol/m3
%Ksp1 = 1.8e-11;& && &% Ksp& && && & mol2/L2
Ksp2 = 3.5e-8;& && & % Ksp& && && && &mol2/L2
k11 = 8.4e3;& && && &% L/(mol*s) k11 = 8.4*10^3& && && && & L/mol*s
k12 = 2e-4;& && && & % s^(-1)& &&&k12 = 2.0*10^(-4)& && && &&&s^(-1)
k21 = 6e9;& && && &&&% L/(mol*s) k21 = 6*10^(9)& && && && & L/mol*s
k22 = 1.2e6;& && && &% s^(-1)& &&&k22 = 1.2*10^(6)& && && && &s^(-1)
k31 = 1.4e11;& && &&&% L/(mol*s) k31 = 1.4*10^(11)& && && & L/mol*s
k32 = 1.3e-3;& && &&&% mol/(L*s) k32 = 1.3*10^(-3)& && && & mol/L*s
k41 = 2.4e-2;& && &&&% s^(-1)& &&&k41 = 2.4*10^(-2)& && && & s^(-1)
k42 = 5.7e2;& && && &% L/(mol*s) k42 = 5.7*10^4& && && && & L/mol*s
k52 = 9e-3;& && && & % s^(-1) k52 = 9e-3& && && && && && && &s^(-1)
k51 = 1.88e6;& && && &%k52/Ksp2;& && &% L/(mol*s)&&
%D_CO2_l&&= 1.91e-7;&&%diffusion coefficeint of CO2(l) dm2/s
D_OH& &&&= 5.27e-7;&&%diffusion coefficeint of OH- dm2/s
D_H& && &= 9.31e-7;&&%diffusion coefficeint of H+ dm2/s
%D_HCO3& &= 1.19e-7;&&%diffusion coefficeint of HCO3 dm2/s
%D_CO3& & = 9.23e-8;&&%diffusion coefficeint of CO3 dm2/s
D_Mg& &&&=1.06e-9; %7.06e-8;&&%diffusion coefficeint of Mg2+ dm2/s
%zk_CO2_l&&= 0;
zk_OH& &&&= -1;
zk_H& && &= 1;
%zk_HCO3& &= -1;
%zk_CO3& & = -2;
zk_Mg& &&&= 2;
%C_CO2_l&&= 0;& && & % concentration on the interface i
C_OH& &&&= 3.3e-4;&&
C_H& && &= 3.03e-11;
%C_HCO3& &= 0;& && &
%C_CO3& & = 0;& && &
C_Mg& &&&= 1.65e-4; % concentration on the interface i
%d_CO2_l& &= C_CO2_l-x(1); % Ci-Cb
d_OH& && &= C_OH-x(3);
d_H& && & = C_H-x(4);
%d_HCO3& & = C_HCO3-x(5);
%d_CO3& &&&= C_CO3-x(6);
d_Mg& && &= C_Mg-x(7); % Ci-Cb
% d_CO2_l& &= x(1)-C_CO2_l; % Cb-Ci
% d_OH& && &= x(3)-C_OH;
% d_H& && & = x(4)-C_H;
% d_HCO3& & = x(5)-C_HCO3;
% d_CO3& &&&= x(6)-C_CO3;
% d_Mg& && &= x(7)-C_Mg; % Cb-Ci
%A_CO2_l& &= (C_CO2_l+x(1))/2; % (Ci+Cb)/2
A_OH& && &= (C_OH+x(3))/2;
A_H& && & = (C_H+x(4))/2;
%A_HCO3& & = (C_HCO3+x(5))/2;
%A_CO3& &&&= (C_CO3+x(6))/2;
A_Mg& && &= (C_Mg+x(7))/2;&&% (Ci+Cb)/2
Ntot1 = 8.17e10;& && &%total number of Mg(OH)2 100g = 8.17e10&&
& && && && && && && &%Ntot1 = 100/(3.14*d^3*2340/6)
UL = 1.01e-1; % dynamic viscosity 1.01e-3 Pas = kg*s/m = 1.01e-1 g*s/dm
VL = 1.01e-4; % kinematic viscosity 1.01e-6 m2/s = 1.01e-4 dm2/s
eT = 14028; % total energy dissipation w/kg = m2/s3*100 = dm2/s3
& && && && && && &
dp = (6*x(8)/(Ntot1*dd_MgOH*3.14))^(1/3); %diameter of solid dm
a = (3.3*x(8)^(2/3)*(3.14*Ntot1)^(1/3)/dd_MgOH^(2/3)/vl);%specific area
% mass transfer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%& &
Re = eT^(1/3)*dp^(4/3)/VL;
%Sc_CO2_l&&= UL/(ddslurry*D_CO2_l);
Sc_OH& &&&= UL/(ddslurry*D_OH);
Sc_H& && &= UL/(ddslurry*D_H);
%Sc_HCO3& &= UL/(ddslurry*D_HCO3);
%Sc_CO3& & = UL/(ddslurry*D_CO3);
Sc_Mg& &&&= UL/(ddslurry*D_Mg);
%Sh_CO2_l&&= (2^5.8+(0.61*(Re^0.58)*(Sc_CO2_l^(1/3)))^5.8)^(1/5.8);
Sh_OH& &&&= (2^5.8+(0.61*(Re^0.58)*(Sc_OH^(1/3)))^5.8)^(1/5.8);
Sh_H& && &= (2^5.8+(0.61*(Re^0.58)*(Sc_H^(1/3)))^5.8)^(1/5.8);
%Sh_HCO3& &= (2^5.8+(0.61*(Re^0.58)*(Sc_HCO3^(1/3)))^5.8)^(1/5.8);
%Sh_CO3& & = (2^5.8+(0.61*(Re^0.58)*(Sc_CO3^(1/3)))^5.8)^(1/5.8);
Sh_Mg& &&&= (2^5.8+(0.61*(Re^0.58)*(Sc_Mg^(1/3)))^5.8)^(1/5.8);
% Parameter fitting!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%Ks_CO2_l&&= Sh_CO2_l*D_CO2_l/
Ks_OH& &&&= Sh_OH*D_OH/dp*X_1;
Ks_H& && &= Sh_H*D_H/dp*X_2;
%Ks_HCO3& &= Sh_HCO3*D_HCO3/
%Ks_CO3& & = Sh_CO3*D_CO3/
Ks_Mg& &&&= Sh_Mg*D_Mg/dp*X_3;
% K = (zk_CO2_l*Ks_CO2_l*d_CO2_l+zk_OH*Ks_OH*d_OH+zk_H*Ks_H*d_H+zk_HCO3*Ks_HCO3*d_HCO3+zk_CO3*Ks_CO3*d_CO3+zk_Mg*Ks_Mg*d_Mg)/...
%& &&&(zk_CO2_l^2*Ks_CO2_l*A_CO2_l+zk_OH^2*Ks_OH*A_OH+zk_H^2*Ks_H*A_H+zk_HCO3^2*Ks_HCO3*A_HCO3+zk_CO3^2*Ks_CO3*A_CO3+zk_Mg^2*Ks_Mg*A_Mg);
K = (zk_OH*Ks_OH*d_OH+zk_H*Ks_H*d_H+zk_Mg*Ks_Mg*d_Mg)/...
& & (zk_OH^2*Ks_OH*A_OH+zk_H^2*Ks_H*A_H+zk_Mg^2*Ks_Mg*A_Mg);
%NCO2_l&&= Ks_CO2_l*(d_CO2_l-A_CO2_l*zk_CO2_l*K);
NOH& &&&= Ks_OH*(d_OH-A_OH*zk_OH*K);
NH& && &= Ks_H*(d_H-A_H*zk_H*K);
%NHCO3& &= Ks_HCO3*(d_HCO3-A_HCO3*zk_HCO3*K);
%NCO3& & = Ks_CO3*(d_CO3-A_CO3*zk_CO3*K);
NMg& &&&= Ks_Mg*(d_Mg-A_Mg*zk_Mg*K);
% K = (zk_CO2_l*D_CO2_l*d_CO2_l+zk_OH*D_OH*d_OH+zk_H*D_H*d_H+zk_HCO3*D_HCO3*d_HCO3+zk_CO3*D_CO3*d_CO3+zk_Mg*D_Mg*d_Mg)/...
%& &&&(zk_CO2_l^2*D_CO2_l*A_CO2_l+zk_OH^2*D_OH*A_OH+zk_H^2*D_H*A_H+zk_HCO3^2*D_HCO3*A_HCO3+zk_CO3^2*D_CO3*A_CO3+zk_Mg^2*D_Mg*A_Mg);
% NCO2_l&&= Ks_CO2_l*d_CO2_l-Ks_CO2_l*A_CO2_l*zk_CO2_l*K;
% NOH& &&&= Ks_OH*d_OH-Ks_OH*A_OH*zk_OH*K;
% NH& && &= Ks_H*d_H-Ks_H*A_H*zk_H*K;
% NHCO3& &= Ks_HCO3*d_HCO3-Ks_HCO3*A_HCO3*zk_HCO3*K;
% NCO3& & = Ks_CO3*d_CO3-Ks_CO3*A_CO3*zk_CO3*K;
% NMg& &&&= Ks_Mg*d_Mg-Ks_Mg*A_Mg*zk_Mg*K;
% reaction kinetics%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1 = k11*x(1)*x(3)-k12*x(5);
r2 = k21*x(5)*x(3)-k22*x(6);
r3 = k31*x(3)*x(4)-k32;
r4 = k41*x(1)-k42*x(5)*x(4);
r5 = k51*x(7)*x(6)-k52;
% Crystallization kinetics%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameter fitting!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% kn = 10^2.6;&&% kn = 10^2.6;%mol^(-n)*dm^(3n-3)*s^-1
% kg = 8.06e-2; % kg = 8.06e-2;%mol^(-g)*dm^(3g+1)*s^-1
% n = 6.5;& && &% n = 4.2;
% g = 6.5;& && &&&% g = 2;
S = sqrt((x(6)+x(9))*(x(7)+x(9)))-sqrt(Ksp2);
B = kn*(S^n);
G = kg*(S^g);
%current = 2*NMg+NH-NOH
dx(1) = kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg)/vl-r1-r4;%+NCO2_l*a;& &
dx(2) = (Qg*(CO2in-x(2))-kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg))/
dx(3) = -r1-r2-r3+NOH*a;
dx(4) = -r3+r4+NH*a;
dx(5) = r1-r2+r4;%+NHCO3*a;
dx(6) = r2-r5;%+NCO3*a;
dx(7) = NMg*a-r5;
dx(8) = -NMg*a*vl*M_MgOH;
dx(9) = r5-pi/6*(1e-21)*dd_mol_MgCO3*B-G/2*pi*dd_mol_MgCO3*x(12); %MgCO3(l) mol/L
dx(10) = B; % m0& &
dx(11) = G*x(10);%m1
dx(12) = 2*G*x(11);%m2
dx(13) = 3*G*x(12);%m3
哦 不好意思 一共有十三个方程组&&但是我只有其中的两组实验数据用于拟合 也就是x(:,7)+x(:,9) 和 -log10(x(:,4)&&只要拟合上这两组，数学模型就可以了&&我在Objfunc 函数中也是只写了这两组数据&&具体您说的*把方程组以完整的文本形式给出* 我不太明白什么意思 我的代码就是从matlab里贴过来了 您可以直接贴到matlab里运行 请问您能说的详细点吗 非常感谢& &另外 我把您先前写的代码在MatlabR2014a里运行 提示lsqnonlin调用 错误 我不太明白 是不是新版本的调用格式有所区别 再次感谢您宝贵的时间
按照我在二楼说的格式，你给出具体方程和数据。因为你的代码我看不懂。
文本形式是指，方程以MATLAB能识别的形式给出，我好直接复制，因为我对你的问题不了解，在调试程序的时候，我来输入方程的话，很可能输错。
我之前写的代码，都是m文件的格式，不要在command window中运行。新建一个m文件，粘贴所有代码，运行即可。
高版本的MATLAB一般都能兼容低版本的。
哦 明白了 我把这个三段代码分别放到 三个m文件里进行调用 所以算不出来， 现在可以运行了 但是每次求出的都是local minimum possible 而且 我给不同的初值 拟合出来的数值都在初值附近 并不是最优解。
关于您说的文本形式，您的意思是说我的第三段ode代码（动力学方程） 您看不懂？你是需要原始的数学模型吗？ 还是第一段lsqnonlin代码不够清楚？ 因为我这个KineticsEqs 包含传质，电荷守恒，吸收，溶解，结晶等模型，很难用简单的y和k表示出来，所以我给你写一个直接能复制的代码，您复制到m文件里直接运行即可（KineticsEqs的参数部分 我用一行%标出 我现在的问题是 如果用nlparci 就会出错 不用nlparci就会出现局部最优解 请您运行下试试 非常感谢）：
function fitting_Wenli_Zhao
format long
% initial parameter identification
% kn = 10^2.6; %Crystallization kinetics
% kg = 8.06e-2;
% n = 6.5 ;
% g = 6.5 ;
k0 = ; %Concentration and PH
%Experimental data
% Concentration and PH
size(data);
%initial values for KineticsEqs
% Boundary of fitting parameter
% Call the optimizer:
options = optimoptions(@lsqnonlin,'TolX',1e-6,'TolFun',1e-6,'MaxFunEvals',500);
& &= lsqnonlin(@ObjFunc,k0,lb,ub,options,tspan,x0,yexp);
%ci = nlparci(teta,residual,jacobian);???????
%Optimized solution
= ode15s(@KineticsEqs,ts,x0,,k);
plot(tspan,yexp(:,1),'o',ts,ys(:,7)+ys(:,9)) %plot the data vs solution
title('Parameter fitting of Mg2+ in the bulk,mol/L');
PH=-log10(ys(:,4));
plot(tspan,yexp(:,2),'o',ts,PH) %plot the data vs solution
title('Parameter fitting of PH');
fprintf('\n\nlsqnonlin():\n')
fprintf('\tkn = %.16f \n',k(1))
fprintf('\tkg = %.16f \n',k(2))
fprintf('\tn&&= %.16f \n',k(3))
fprintf('\tg&&= %.16f \n',k(4))
fprintf('\tf&&= %.16f \n',k(5))
%Call 'Objfun' function:
function lsq = ObjFunc(k,tspan,x0,yexp)
= ode15s(@KineticsEqs,tspan,x0,,k);
%Concentration of Mg2+
Xsim1 = Xsim(:,7)+Xsim(:,9);
Xsim2 = -log10(Xsim(:,4));
ycal(:,1) = Xsim1;
ycal(:,2) = Xsim2;
size(ycal);
size(ycal);
% the reaction of Calc hydroxide and Carbon dioxide
function dxdt=KineticsEqs(t,x,k)
dxdt = zeros(13,1);
% x(1) is the concentration of liquid CO2(l)
% initial value = 0 mol/L
% x(2) is the concentration of gas CO2(g)
% initial value = 0 mol/L
% x(3) is the concentration of Hydroxy&&in bulk
% initial value = 3.3e-4 mol/L
% x(4) is the concentration of Hydrogen&&in bulk
% initial value = 3.03e-11 mol/L
% x(5) is the concentration of Bicarbonate&&
% initial value = 0 mol/L
% x(6) is the concentration of carbonate&&
% initial value = 0 mol/L
% x(7) is the concentration of Magnesium
% initial value = 1.65e-4&&mol/L
% x(8) is the mass of Mg(OH)2
% initial value = 100&&g
% x(9) is the concentraiton of liquid MgCO3
% initial value = 0&&mol/L
% x(10) is o monent which is number density of crystal n(L)
% initial value = 0
% x(11) is 1 monent which is number density of crystal n(L)
% initial value = 0
% x(12) is 2 monent which is number density of crystal n(L)
% initial value = 0
% x(13) is 3 monent which is number density of crystal n(L)
% initial value = 0
Qg = 1.67e-2;& && &&&% L/s& && && & 1L/min
N&&= 650;&&
vl = 2.826;& && && & % vl is the dispersion volum& && && && && && &L
vg = 2.047*10^(-5)*N^(0.579)*Qg^0.63*&&
& && && && && && && &% vg is the gas volum in the gas-liquid mixture m3
& && && &&&%agitating speed rpm
P = 1;& && && && && &% P is the atmospheric pressure& && && && && & atm
pi = 3.14;
kla = 0.0304;& && &&&% kla is the mass transfer coefficient of CO2&&s^(-1)
& && && && && && && &%(9.6218e-5)*N^0.4316*Qg^0.3840*c^0.1117*E;&&
kH = 0.035;& && && & % kh is the Henry's constant& && && && && &mol/(L*atm)
M_CO2 = 44;& && && & % M_CO2 is the co2 molar mass& && && && &&&g/mol
M_MgOH = 58;& && && &% molar mass of Mg(OH)2& && && && && &g/mol
M_MgCO3 = 84;& && &&&% molar mass of M_MgCO3& & g/mol
dd_CO2 = 1.8;& && &&&% dd_CO2 is the density of CO2
dd_MgOH = 2.34e3;& & % dd_Mg(OH)2 is the density of Mg(OH)2& && && & g/L
dd_MgCO3 = 2.96e3;& &% dd_Mg(OH)2 is the density of MgCO3& && && & g/L
dd_mol_MgCO3 = dd_MgCO3/M_MgCO3; % mol/L
ddslurry = 1e3;& && &% density of slurry g/L& && && && && && && &% ddCO2 = 1960; ddCO2 is the density of CO2
c = 1;& && && && && &% volume fraction of CO2 in the inlet gas& && && && && && && &
CO2in = dd_CO2/M_CO2*c;% CO2in is the concentration of CO2 entering the reactor mol/L
& && && && && && && &% CO2in = ddCO2/MCO2*0.15 =.05 =6.75&&mol/m3
%Ksp1 = 1.8e-11;& && &% Ksp& && && & mol2/L2
Ksp2 = 3.5e-8;& && & % Ksp& && && && &mol2/L2
k11 = 8.4e3;& && && &% L/(mol*s) k11 = 8.4*10^3& && && && & L/mol*s
k12 = 2e-4;& && && & % s^(-1)& &&&k12 = 2.0*10^(-4)& && && &&&s^(-1)
k21 = 6e9;& && && &&&% L/(mol*s) k21 = 6*10^(9)& && && && & L/mol*s
k22 = 1.2e6;& && && &% s^(-1)& &&&k22 = 1.2*10^(6)& && && && &s^(-1)
k31 = 1.4e11;& && &&&% L/(mol*s) k31 = 1.4*10^(11)& && && & L/mol*s
k32 = 1.3e-3;& && &&&% mol/(L*s) k32 = 1.3*10^(-3)& && && & mol/L*s
k41 = 2.4e-2;& && &&&% s^(-1)& &&&k41 = 2.4*10^(-2)& && && & s^(-1)
k42 = 5.7e2;& && && &% L/(mol*s) k42 = 5.7*10^4& && && && & L/mol*s
k52 = 9e-3;& && && & % s^(-1) k52 = 9e-3& && && && && && && &s^(-1)
k51 = 1.88e6;& && && &%k52/Ksp2;& && &% L/(mol*s)&&
%D_CO2_l&&= 1.91e-7;&&%diffusion coefficeint of CO2(l) dm2/s
D_OH& &&&= 5.27e-7;&&%diffusion coefficeint of OH- dm2/s
D_H& && &= 9.31e-7;&&%diffusion coefficeint of H+ dm2/s
%D_HCO3& &= 1.19e-7;&&%diffusion coefficeint of HCO3 dm2/s
%D_CO3& & = 9.23e-8;&&%diffusion coefficeint of CO3 dm2/s
D_Mg& &&&=1.06e-9; %7.06e-8;&&%diffusion coefficeint of Mg2+ dm2/s
%zk_CO2_l&&= 0;
zk_OH& &&&= -1;
zk_H& && &= 1;
%zk_HCO3& &= -1;
%zk_CO3& & = -2;
zk_Mg& &&&= 2;
%C_CO2_l&&= 0;& && & % concentration on the interface i
C_OH& &&&= 3.3e-4;&&
C_H& && &= 3.03e-11;
%C_HCO3& &= 0;& && &
%C_CO3& & = 0;& && &
C_Mg& &&&= 1.65e-4; % concentration on the interface i
%d_CO2_l& &= C_CO2_l-x(1); % Ci-Cb
d_OH& && &= C_OH-x(3);
d_H& && & = C_H-x(4);
%d_HCO3& & = C_HCO3-x(5);
%d_CO3& &&&= C_CO3-x(6);
d_Mg& && &= C_Mg-x(7); % Ci-Cb
% d_CO2_l& &= x(1)-C_CO2_l; % Cb-Ci
% d_OH& && &= x(3)-C_OH;
% d_H& && & = x(4)-C_H;
% d_HCO3& & = x(5)-C_HCO3;
% d_CO3& &&&= x(6)-C_CO3;
% d_Mg& && &= x(7)-C_Mg; % Cb-Ci
%A_CO2_l& &= (C_CO2_l+x(1))/2; % (Ci+Cb)/2
A_OH& && &= (C_OH+x(3))/2;
A_H& && & = (C_H+x(4))/2;
%A_HCO3& & = (C_HCO3+x(5))/2;
%A_CO3& &&&= (C_CO3+x(6))/2;
A_Mg& && &= (C_Mg+x(7))/2;&&% (Ci+Cb)/2
Ntot1 = 8.17e10;& && &%total number of Mg(OH)2 100g = 8.17e10&&
& && && && && && && &%Ntot1 = 100/(3.14*d^3*2340/6)
UL = 1.01e-1; % dynamic viscosity 1.01e-3 Pas = kg*s/m = 1.01e-1 g*s/dm
VL = 1.01e-4; % kinematic viscosity 1.01e-6 m2/s = 1.01e-4 dm2/s
eT = 14028; % total energy dissipation w/kg = m2/s3*100 = dm2/s3
& && && && && && &
dp = (6*x(8)/(Ntot1*dd_MgOH*3.14))^(1/3); %diameter of solid dm
a = (3.3*x(8)^(2/3)*(3.14*Ntot1)^(1/3)/dd_MgOH^(2/3)/vl);%specific area
% mass transfer& && &
Re = eT^(1/3)*dp^(4/3)/VL;
%Sc_CO2_l&&= UL/(ddslurry*D_CO2_l);
Sc_OH& &&&= UL/(ddslurry*D_OH);
Sc_H& && &= UL/(ddslurry*D_H);
%Sc_HCO3& &= UL/(ddslurry*D_HCO3);
%Sc_CO3& & = UL/(ddslurry*D_CO3);
Sc_Mg& &&&= UL/(ddslurry*D_Mg);
%Sh_CO2_l&&= (2^5.8+(0.61*(Re^0.58)*(Sc_CO2_l^(1/3)))^5.8)^(1/5.8);
Sh_OH& &&&= (2^5.8+(0.61*(Re^0.58)*(Sc_OH^(1/3)))^5.8)^(1/5.8);
Sh_H& && &= (2^5.8+(0.61*(Re^0.58)*(Sc_H^(1/3)))^5.8)^(1/5.8);
%Sh_HCO3& &= (2^5.8+(0.61*(Re^0.58)*(Sc_HCO3^(1/3)))^5.8)^(1/5.8);
%Sh_CO3& & = (2^5.8+(0.61*(Re^0.58)*(Sc_CO3^(1/3)))^5.8)^(1/5.8);
Sh_Mg& &&&= (2^5.8+(0.61*(Re^0.58)*(Sc_Mg^(1/3)))^5.8)^(1/5.8);
% Parameter fitting!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%Ks_CO2_l&&= Sh_CO2_l*D_CO2_l/
Ks_OH& &&&= Sh_OH*D_OH/
Ks_H& && &= Sh_H*D_H/
%Ks_HCO3& &= Sh_HCO3*D_HCO3/
%Ks_CO3& & = Sh_CO3*D_CO3/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%k(5)
Ks_Mg& &&&= Sh_Mg*D_Mg/dp*k(5);
% K = (zk_CO2_l*Ks_CO2_l*d_CO2_l+zk_OH*Ks_OH*d_OH+zk_H*Ks_H*d_H+zk_HCO3*Ks_HCO3*d_HCO3+zk_CO3*Ks_CO3*d_CO3+zk_Mg*Ks_Mg*d_Mg)/...
%& &&&(zk_CO2_l^2*Ks_CO2_l*A_CO2_l+zk_OH^2*Ks_OH*A_OH+zk_H^2*Ks_H*A_H+zk_HCO3^2*Ks_HCO3*A_HCO3+zk_CO3^2*Ks_CO3*A_CO3+zk_Mg^2*Ks_Mg*A_Mg);
K = (zk_OH*Ks_OH*d_OH+zk_H*Ks_H*d_H+zk_Mg*Ks_Mg*d_Mg)/...
& & (zk_OH^2*Ks_OH*A_OH+zk_H^2*Ks_H*A_H+zk_Mg^2*Ks_Mg*A_Mg);
%NCO2_l&&= Ks_CO2_l*(d_CO2_l-A_CO2_l*zk_CO2_l*K);
NOH& &&&= Ks_OH*(d_OH-A_OH*zk_OH*K);
NH& && &= Ks_H*(d_H-A_H*zk_H*K);
%NHCO3& &= Ks_HCO3*(d_HCO3-A_HCO3*zk_HCO3*K);
%NCO3& & = Ks_CO3*(d_CO3-A_CO3*zk_CO3*K);
NMg& &&&= Ks_Mg*(d_Mg-A_Mg*zk_Mg*K);
% K = (zk_CO2_l*D_CO2_l*d_CO2_l+zk_OH*D_OH*d_OH+zk_H*D_H*d_H+zk_HCO3*D_HCO3*d_HCO3+zk_CO3*D_CO3*d_CO3+zk_Mg*D_Mg*d_Mg)/...
%& &&&(zk_CO2_l^2*D_CO2_l*A_CO2_l+zk_OH^2*D_OH*A_OH+zk_H^2*D_H*A_H+zk_HCO3^2*D_HCO3*A_HCO3+zk_CO3^2*D_CO3*A_CO3+zk_Mg^2*D_Mg*A_Mg);
% NCO2_l&&= Ks_CO2_l*d_CO2_l-Ks_CO2_l*A_CO2_l*zk_CO2_l*K;
% NOH& &&&= Ks_OH*d_OH-Ks_OH*A_OH*zk_OH*K;
% NH& && &= Ks_H*d_H-Ks_H*A_H*zk_H*K;
% NHCO3& &= Ks_HCO3*d_HCO3-Ks_HCO3*A_HCO3*zk_HCO3*K;
% NCO3& & = Ks_CO3*d_CO3-Ks_CO3*A_CO3*zk_CO3*K;
% NMg& &&&= Ks_Mg*d_Mg-Ks_Mg*A_Mg*zk_Mg*K;
% reaction kinetics
r1 = k11*x(1)*x(3)-k12*x(5);
r2 = k21*x(5)*x(3)-k22*x(6);
r3 = k31*x(3)*x(4)-k32;
r4 = k41*x(1)-k42*x(5)*x(4);
r5 = k51*x(7)*x(6)-k52;
% Crystallization kinetics
% Parameter fitting!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% kn = 10^2.6;&&% kn = 10^2.6;%mol^(-n)*dm^(3n-3)*s^-1
% kg = 8.06e-2; % kg = 8.06e-2;%mol^(-g)*dm^(3g+1)*s^-1
% n = 6.5;& && &% n = 4.2;
% g = 6.5;& && &&&% g = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k(1)-k4
S = sqrt((x(6)+x(9))*(x(7)+x(9)))-sqrt(Ksp2);
B = k(1)*(S^k(3));
G = k(2)*(S^k(4));
%current = 2*NMg+NH-NOH
dxdt(1) = kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg)/vl-r1-r4;%+NCO2_l*a;& &
dxdt(2) = (Qg*(CO2in-x(2))-kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg))/
dxdt(3) = -r1-r2-r3+NOH*a;
dxdt(4) = -r3+r4+NH*a;
dxdt(5) = r1-r2+r4;%+NHCO3*a;
dxdt(6) = r2-r5;%+NCO3*a;
dxdt(7) = NMg*a-r5;
dxdt(8) = -NMg*a*vl*M_MgOH;
dxdt(9) = r5-pi/6*(1e-21)*dd_mol_MgCO3*B-G/2*pi*dd_mol_MgCO3*x(12); %MgCO3(l) mol/L
dxdt(10) = B; % m0& &
dxdt(11) = G*x(10);%m1
dxdt(12) = 2*G*x(11);%m2
dxdt(13) = 3*G*x(12);%m3
你的代码我基本看懂了。计算结果local minimum possible，是正常的，因为本身lsqnonlin就是局部最优的算法，这种显示表明计算的步骤能跑完。ci = nlparci(teta,residual,jacobian)这是用算拟合参数的置信区间的，这一块可以先不管，先把拟合结果调试好了再说。
在确保程序中参数都输入正确、程序可以跑通的情况下，主要精力就放在调初值上面，如上面所说lsqnonlin的局限性，对于参数个数较多的情况，用lsqnonlin获得很好的结果还是比较困难的。
建议你在程序中加一个计算决定系数的代码，具体参见：
http://emuch.net/bbs/viewthread.php?tid=7603645，其中的“R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);”，先把R2调到0.9以上再说。
非常感谢！我刚从外面开会回来所以回复的比较晚 不要意思。 我这个拟合效果不好 也有一部分是模型的原因 因为这个是我自己建的数学模型 所以在模型的准确性上还存在一些问题。 另外 我现在遇到的难题是 拟合可以跑通并得到参数值 但是如果把这个参数带到原来的ode去解 解出来的结果和拟合是解出来的结果不一样。 我自己对比了两个代码，因为是复制过去的不可能存在不同。请问这是拟合中的常见问题么？谢谢
您好 我上述模型问题解决了 但是对于lsqnonlin&&fitting
我无论怎么变化初值 得到的结果 都还是给的初值 并没有实现优化的过程 我不知道code哪出问题了 您能帮我运行下试试么 非常感谢
建议你联系一下 dingd，请他用1stopt帮你跑一下。
研究生必备与500万研究生在线互动！
扫描下载送金币