求会matlab的化工大佬帮我看一下我动力学模型的问题
function nihe_0316% 微分方程的参数估计
%
clear all;
clc;
format long;
global t0;
% k10,E,Ka 对应k(1)、k(2)、k(3)
k0 =; % 参数初值，重要！，可将结果反复代入迭代计算！！
x0 = 0.63713;%初始状态，没有就取第一行数据
lb = ; % 参数下限,根据拟合结果，适当调整
ub = [+inf +inf +inf]; % 参数上限
% 实验数据:
expdata= ...
60 0.79642
120 0.87377
240 0.90782
360 0.91735];
t0=expdata(:,1);
yexp =expdata(:,2); % yexp: 对应实验数据
% ------------------------------------------------------------------
% 使用函数lsqnonlin()进行参数估计
k=lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk10 = %.8f\n',k(1))
fprintf('\tE = %.8f\n',k(2))
fprintf('\tKa = %.8f\n',k(3))
% ----------------------------------------------------
tspan = ;
= ode45(@KineticEqs,tspan,x0,[],k);
x=spline(t,xx,t0); %插值，使理论值和实验值个数一致
% 计算相关系数
A=corrcoef(yexp,x(:,1));
fprintf('\nXb的相关系数R为:\n')
fprintf('\tR = %.8f\n',A(1,2))
%
% 计算决定系数
r(1)=1-sum((yexp-x(:,1)).^2)./sum((yexp-mean(yexp)).^2);
% ********************************决定系数
fprintf('\nXb的决定系数R^2为:\n')
fprintf('\tR^2 = %.8f\n',r(1))
% 计算均方根误差
rmse(1)=sqrt(sum((yexp-x(:,1)).^2)/length(t0));
fprintf('\nXb的均方根误差RMSE为:\n')
fprintf('\tRMSE = %.8f\n',rmse(1))
%
% 计算残差平方和 SSE（Sum of Squares for Error）
ssr(1)=sum((yexp(:,1)-x(:,1)).^2);
fprintf('\nXb的残差平方和为:\n')
fprintf('\tSSE = %.8f\n',ssr(1))
% ----------------------------------------------------
figure(1)
plot(t,xx(:,1),'b-',t0,yexp,'ro'),legend('理论值','实验值','Location','best')
xlabel('时间');ylabel('实验值N');
% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
tspan = ;
= ode45(@KineticEqs,tspan,x0,[],k);
x=spline(t,xx,t0);
f = x - yexp; % 理论值与实验值的差值，残差
end
% ------------------------------------------------------------------
function dXbdt = KineticEqs(t,Xb,k) % 微分方程, k10,E,Ka 对应k(1)、k(2)、k(3)
T=358.15;
m=3;
N0=0.1467;
Cb0=0.5564;
Ca=Cb0*(3-Xb);
Cb=Cb0*(1-Xb);
Cc=Cb0*Xb;
Keq=exp(4873/T-13.813);
Kb=5.2112*k(3);
Kc=2.2695*k(3);
dXbdt = (m*k(1)*exp(-k(2)/8.314*T)*(Ca*Cb-Cc/Keq))/(N0*(1+k(3)*Ca+Kb*Cb+Kc*Cc)^2);
end
% ------------------------------------------------------------------
end
结果算出来的X的计算值不随时间变化，理论上是一个dx/dt=f（x）的方程，我看书上dc/dt=f（c）的方程用RK法可以计算加拟合，不知道为什么我这个就不行，请大家帮我看一下
学习学习，谢谢分享 学习学习，谢谢分享
页:
[1]