["

Thread Content

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 = ; % 参数上限 % 实验数据: expdata= ... ; 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法可以计算加拟合，不知道为什么我这个就不行，请大家帮我看一下