常微分方程拟合matlab,常微分方程组参数拟合问题，跪求啊，两年了，各位大神帮帮忙啊...

各位论坛里的大神帮帮忙啊，让我早日出坑吧

因为word里的公式粘贴不过来，只能把问题添加到附件了，各位大神不要闲麻烦啊，点开看看哈

自己编的程序如下;

function Kinetics4

clear all

clc

k0 = [2.159*10^8 4.734*10^6 174686 149892 0.99 0.7 0.5 0.78];%参数初值，其实是不知道的

lb = [-inf -inf -inf -inf -inf -inf -inf -inf];                   % 参数下限

ub = [inf inf inf inf inf inf inf inf];    % 参数上限

u0 = [0,0]; %y初值

a=[0;0.06680583;0.192617534;0.301693824;0.419783943;0.527041818;0.598345407;0.65055108;0.6876854;0;0.0295189820000000;0.0760025420000000;0.143652564000000;0.303435528000000;0.419715167000000;0.525496977000000;0.565069320000000;0.743536298000000;0;0.0411725930000000;0.0608472370000000;0.0986899250000000;0.241132264000000;0.324760943000000;0.466871464000000;0.531688237000000;0.683614442000000];

b=[0;0.2378;0.3288;0.3556;0.39;0.3877;0.3766;0.3580;0.3384;0;0.1373;0.2123;0.2755;0.3561;0.3549;0.3475;0.3486;0.2518;0;0.115879961000000;0.178182567000000;0.219982937000000;0.315058656000000;0.339244190000000;0.348664305000000;0.351963292000000;0.313743618000000];

yexp=[a,b];

% 使用函数lsqnonlin()进行参数估计

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],u0,yexp);

ci = nlparci(k,residual,jacobian);

fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))

fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))

fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))

fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))

fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))

fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))

fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))

fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))

fprintf('  The sum of the squares is: %.1e\n\n',resnorm)

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')

%目标函数

function f =ObjFunc4LNL(k,u0,yexp)

global theta Pt Ac R T Fa0

theta=zeros(5,1);

r0= 1.24*10^-3; % Outer diameter  of the membrane, [m]

ri= 9.4*10^-4;% Inner diameter  of the membrane, [m]

r=(r0-ri)/log(r0/ri);% equivalent diameter ,[m]

Pi=3.14159;

L=0.05;% membrane length,[m]

Ac=Pi*r*L;% membrane area [m^2]

R=8.314; % [J/K.mol]

Pt=101.325;

Ftot=[2.678*10^-2;2.678*10^-2*1.5;2.678*10^-2*2];

Kmax=length(Ftot);

T0=600+273;% Inlet temperature [K]

nmax=9;

X1=zeros(nmax,Kmax);X2=zeros(nmax,Kmax);

for K=1:Kmax

Fa0=Ftot(K); % If Pt is the parameter/Matrix

%Pt=Pt0+(k-1)*1; % If Pt is an axis

%L=Lsize(k); % If L is the parameter/Matrix

%theta(2)=H2Oratio+(k-1)*1; % If m is the parameter/Matrix

theta(2)=3; % F0(H2O)/F0(CH4)

theta(3)=0; % F0(CO)/F0(CH4)

theta(4)=0; % F0(CO2)/F0(CH4)

theta(5)=0; % F0(H2)/F0(CH4)

for n=1:nmax

%delta=delta0+(n-1)*1e-6; % If delta is the parameter/Matrix

%rit=r0+delta; % inner tube radius [m]

%sweep=s0+(n-1); % If sweep is the parameter/Matrix

%T=T0; % If T is a constant

T=T0+(n-1)*50;

ksispan=(0:0.1:1); % precision

[ksi,S1]=ode23s(@KineticEqs,ksispan,u0,[],k); % PBR reactor

X1(n,Kmax)=real(S1(end,1)); % methane finale conversion

X2(n,Kmax)=real(S1(end,2)); % carbon dioxide finale conversion

end

end

f1=X1(n,Kmax)-yexp(:,1);

f2=X2(n,Kmax)-yexp(:,2);

f=[f1;f2];

%T0=T0-273.15;T=T-273.15; % Temperature expressed in Celsius

%Tspan=(T0:1:T);

%Pspan=(Pt0:1:Pt);

%Dspan=(delta0*1e6:1:delta*1e6);

%Sspan=(s0:1:sweep);

%D(:,:)=X3(:,:)-X1(:,:); % Delta = XCH4(MR) - XCH4(PBR)

function diff=KineticEqs(ksi,u,k)

global theta  Pt  Ac Fa0 T R

S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));

P(1)=(1-u(1)-u(2))*S;

P(2)=(theta(2)-u(1)-2*u(2))*S;

P(3)=(theta(3)+u(1))*S;

P(4)=(theta(4)+u(2))*S;

P(5)=(theta(5)+3*(u(1)+u(2))+u(2))*S;

% Differential System

diff1=Ac/Fa0*k(1)*exp(-k(3)/R/T)*P(1)^k(5)*P(2)^k(7);

diff2=Ac/Fa0*k(2)*exp(-k(4)/R/T)*P(1)^k(6)*P(2)^k(8);

diff=[diff1 diff2]';

运行结果：

In nlparci (line 104)

In Kinetics4 (line 15)

警告: 矩阵为奇异工作精度。

k1 = 215900000.0000 ± NaN

k2 = 4734000.0000 ± NaN

k3 = 174686.0000 ± NaN

k4 = 149892.0000 ± NaN

k5 = 1.2381 ± NaN

k6 = 0.5426 ± NaN

k7 = -0.3426 ± NaN

k8 = 0.4455 ± NaN

The sum of the squares is: 2.0e+00

2020-12-6 16:38 上传

点击文件名下载附件

100 KB, 下载次数: 14