Thread Content
funçãonihe_0316 % Estimativa de parâmetros da equação diferencial% limpar tudo; clc; formato longo; global t0; % k10,E,Ka corresponde a k(1), k(2), k(3) k0 =; % valor inicial do parâmetro, importante! , os resultados podem ser substituídos repetidamente em cálculos iterativos! ! x0 = 0,63713; %Estado inicial, caso contrário, pegue a primeira linha de dados lb = ; % Limite inferior do parâmetro, ajuste adequadamente de acordo com os resultados do ajuste ub = ; % Limite superior do parâmetro% Dados experimentais: expdata= ... ; t0=expdados(:,1); yexp =expdata(:,2); % yexp: Correspondente aos dados experimentais% ------------------------------------------------------------------------------- % Use a função lsqnonlin() para estimativa de parâmetros k=lsqnonlin(@ObjFunc4LNL,k0,lb,ub,,x0,yexp); fprintf('\n\nO valor do parâmetro estimado usando a função 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); %Interpolação para tornar o valor teórico consistente com o valor experimental% Calcular o coeficiente de correlação A=corrcoef(yexp,x(:,1)); fprintf('\nO coeficiente de correlação R de% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Coeficiente de determinação fprintf('\nO coeficiente de determinação R^2 de fprintf('\tRMSE = %.8f\n',rmse(1)) % % Calcule a soma residual dos quadrados SSE (soma dos quadrados para erro) ssr(1)=sum((yexp(:,1)-x(:,1)).^2); fprintf('\nA soma residual dos quadrados de Xb é:\n') fprintf('\tSSE = %.8f\n',ssr(1)) % ---------------------------------------------------------------- figure(1) plot(t,xx(:,1),'b-',t0,yexp,'ro'),legend('valor teórico','valor experimental','Localização','melhor') xlabel('tempo');ylabel('valor experimental N'); ObjFunc4LNL(k,x0,yexp) tspan = ; = ode45(@KineticEqs,tspan,x0,,k); * (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 O valor calculado de