MATLABが使える化学工学の専門家の方、この動力学モデルに問題がないか見ていただけませんか

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の計算 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法を用いて計算・フィッティングできると書かれている。なぜ私の場合はうまくいかないのか、皆さんに見てもらいたい