Ich bitte jemanden, der mit MATLAB umgehen kann und aus dem Bereich Chemie ist, mir zu helfen, zu überprüfen, ob mit diesem Dynamikmodell etwas nicht in Ordnung ist

Thread Content

function nihe_0316 % Parameterabschätzung für Differentialgleichungen % clear all; clc; format long; global t0; % k10,E,Ka entsprechen k(1), k(2), k(3) k0 =; % Anfangswerte der Parameter, wichtig! Die Ergebnisse können wiederholt in die iterative Berechnung eingesetzt werden! ! x0 = 0,63713; % Anfangswert – falls nicht vorhanden, werden die Daten der ersten Zeile verwendet.
lb = ; % Untergrenze der Parameter; muss anhand der Ergebnisse der Anpassung angepasst werden.
ub = ; % Obergrenze der Parameter.
% Experimentelle Daten: expdata= … ; t0=expdata(:,1); yexp = expdata(:,2); % yexp: Entspricht den experimentellen Daten.
% ------------------------------------------------------------------
% Nutzung der Funktion lsqnonlin() zur Schätzung der Parameter:
k = lsqnonlin(@ObjFunc4LNL, k0, lb, ub, x0, yexp);
fprintf(“\n\nDie durch die Funktion lsqnonlin() ermittelten Parameterwerte sind:\n”)
fprintf(“\tk10 = %.8f\n”, k(1))
fprintf(“\tE = %.8f\n”, k(2))
fprintf(“\tKa = %.8f\n”, k(3))
% ----------------------------------------------------
tspan = ;
x = ode45(@KineticEqs, tspan, x0, k);
x = spline(t, xx, t0); % Interpolation, damit die Anzahl der theoretischen und experimentellen Werte übereinstimmt.
% Berechnung des Korrelationskoeffizienten:
A = corrcoef(yexp, x(:,1));
fprintf(“\nDer Korrelationskoeffizient von Xb beträgt:\n”)
fprintf(“\tR = %.8f\n”, A(1,2))
% Berechnung des Determinationskoeffizienten:
r(1) = 1 – sum((yexp – x(:,1))^2) / sum((yexp – mean(yexp))^2);
% ******************************** Determinationskoeffizient
fprintf(“\nDer Determinationskoeffizient R^2 von Xb beträgt:\n”)
fprintf(“\tR^2 = %.8f\n”, r(1))
% Berechnung des quadratischen Mittelwerts der Fehler:
rmse(1) = sqrt(sum((yexp – x(:,1))^2) / length(t0));
fprintf(“\nDer quadratische Mittelwert der Fehler RMSE von Xb beträgt:\n”)
fprintf(“\tRMSE = %.8f\n”, rmse(1))
% Berechnung der Summe der Quadraten der Fehler:
ssr(1) = sum((yexp(:,1) – x(:,1))^2);
fprintf(“\nDie Summe der Quadraten der Fehler SSE von Xb beträgt:\n”)
fprintf(“\tSSE = %.8f\n”, ssr(1))
% ----------------------------------------------------
figure(1); plot(t, xx(:,1), ‘b-’, t0, yexp, ‘ro’); legend(‘Theoretische Werte’, ‘Experimentelle Werte’, ‘Location’, ‘best’); xlabel(‘Zeit’); ylabel(‘Experimentelle Werte N’);
% ------------------------------------------------------------------
function f = ObjFunc4LNL(k, x0, yexp)
tspan = ;
x = ode45(@KineticEqs, tspan, x0, k);
x = spline(t, xx, t0);
f = x – yexp; % Differenz zwischen theoretischen und experimentellen Werten, also der Fehler.
end
% ------------------------------------------------------------------
function dXbdt = KineticEqs(t, Xb, k)
% Differentialgleichung; k10, E, Ka entsprechen 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
% ------------------------------------------------------------------
Ende.
Die berechneten Werte für X ändern sich mit der Zeit nicht. Theoretisch handelt es sich dabei um eine Gleichung der Form dx/dt = f(x). In den Büchern wird gesagt, dass man für Gleichungen der Form dc/dt = f(c) die RK-Methode verwenden kann, um Lösungen zu finden. Ich verstehe nicht, warum das bei mir nicht funktioniert. Bitte helft mir dabei!