ode45 求解 diff.equation 并进一步拟合 exp.results

Question

我正在构建一个代码来解决差异。等式：

function dy = KIN1PARM(t,y,k)
%
% version : first order reaction
%   A  --> B
%   dA/dt = -k*A
%   integrated form A = A0*exp(-k*t)
%
dy = -k.*y; 
end  

我希望对这个方程进行数值求解，并将结果（y 作为 t 和 k 的函数）用于最小化实验值以获得参数 k 的最佳值。

function SSE = SSE_minimization_1parm(tspan_inp,val_exp,k_inp,y0_inp)
         f = @(Tt,Ty) KIN1PARM(Tt,Ty,k_inp);  %function to call ode45
         size_limit = length(y0_inp);
         options = odeset('NonNegative',1:size_limit,'RelTol',1e-4,'AbsTol', 1e-4);
         [ts,val_theo] = ode45(f, tspan_inp, y0_inp,options);  %Cexp is the state variable     predicted by the model
         err = val_exp - val_theo;
         SSE = sum(err.^2);   %sum squared-error

绘制实验和计算数据的主要代码是：

% Analyzing first order kinetics
clear all; clc;
figure_title = 'Experimental Data';
label_abscissa = 'Time [s]';
label_ordinatus = 'Concentration [mol/L]';
 %
 abscissa = [  0;
            240;
            480;
            720;
            960;
            1140;
            1380;
            1620;
            1800;
            2040;
            2220;
            2460;
            2700;
            2940];
ordinatus = [  0;
            19.6;
            36.7;
            49.0;
            57.1;
            64.5;
            71.4;
            75.2;
            78.7;
            81.3;
            83.3;
            85.5;
            87.0;
            87.7];
%
 title_string = ['  Time [s]', '  |  ', ' Complex [mol/L] ', ' '];
disp(title_string);
for i=1:length(abscissa)
            report_raw_data{i} = sprintf('%1.3E\t',abscissa(i),ordinatus(i));
            disp([report_raw_data{i}]);
end;
%---------------------/plotting dot data/------------------------------------- 
%
f = figure('Position', [100 100 700 500]);
title(figure_title,'FontName','arial','FontWeight','bold', 'FontSize', 12);
xlabel(label_abscissa, 'FontSize', 12);
ylabel(label_ordinatus, 'FontSize', 12);
%
grid on; hold on;
%
marker_style = { 's'};
%
plot(abscissa,ordinatus, marker_style{1},... 
                                'MarkerFaceColor', 'black',...
                                'MarkerEdgeColor', 'black',...
                                'MarkerSize',4);
%---------------------/Analyzing/----------------------------------------
%
options = optimset('Display','iter','TolFun',1e-4,'TolX',1e-4);
%
        CPUtime0 = cputime;
        Time_M = abscissa;
        Concentration_M = ordinatus;
        tspan = Time_M;
        y0 = 0;
        k0 = rand(1);
[k, fval, exitflag, output] = fminsearch(@(k)      SSE_minimization_1parm(tspan,Concentration_M,k,y0),k0,options);
        CPUtimex = cputime;
        CPUtime_delay = CPUtimex - CPUtime0;
%        
%---------------------/plotting calculated data/-------------------------------------
%
xupperlimit = Time_M(length(Time_M));
xval = ([0:1:xupperlimit])';
%
yvector = data4plot_1parm(xval,k,y0);
plot(xval,yvector, 'r');
hold on;
%---------------------/printing calculated data/-------------------------------------
%
disp('RESULTS:');
disp(['CPU time:    ',sprintf('%0.5f\t',CPUtime_delay),' sec']);
disp(['k:       ',sprintf('%1.3E\t',k')]);
disp(['fval:        ',sprintf('%1.3E\t',fval)]);
disp(['exitflag:   ',sprintf('%1.3E\t',exitflag)]);
disp(output);
disp(['Output:      ',output.message]);

对应的函数，它使用优化的参数k来产生计算出的y = f(t)数据：

function val = data4plot_1parm(tspan_inp,k_inp,y0_inp)
         f = @(Tt,Ty) KIN1PARM(Tt,Ty,k_inp);  
         size_limit = length(y0_inp);
         options = odeset('NonNegative',1:size_limit,'RelTol',1e-4,'AbsTol',1e-4);
        [ts,val_theo] = ode45(f, tspan_inp, y0_inp, options); 

代码运行优化周期时总是给出不同的参数 k 值，这与使用 ln(y) vs t 计算的值不同（对于该系列的 exp.数据，应该在 7.0e-4 左右）。

查看 ode 求解器 (SSE_minimization_1parm => val_theo) 的结果，我发现 ode 函数给了我一个零向量。

有人可以帮我弄清楚 ode 求解器的情况吗？

非常感谢！

---

Answer 1

所以我现在能得到的最好的东西来了。以我的方式，我将纵坐标值作为时间，将横坐标值作为您尝试建模的测量量。此外，您似乎为求解器设置了很多选项，我都省略了。首先是您提出的使用ode45() 的解决方案，但使用非零y0 = 100，我只是通过查看数据“猜测”（在半对数图中）。

function main 

abscissa = [0; 
            240;
            480;
            720;
            960;
            1140;
            1380;
            1620;
            1800;
            2040;
            2220;
            2460;
            2700;
            2940];

ordinatus = [  0;
            19.6;
            36.7;
            49.0;
            57.1;
            64.5;
            71.4;
            75.2;
            78.7;
            81.3;
            83.3;
            85.5;
            87.0;
            87.7];

tspan = [min(ordinatus), max(ordinatus)]; % // assuming ordinatus is time

y0 = 100; % // <---- Probably the most important parameter to guess
k0 = -0.1; % // <--- second most important parameter to guess (negative for growth)

        k_opt = fminsearch(@minimize, k0) % // optimization only over k
        % nested minimization function
        function e = minimize(k)
            sol = ode45(@KIN1PARM, tspan, y0, [], k);
            y_hat = deval(sol, ordinatus); % // evaluate solution at given times
            e = sum((y_hat' - abscissa).^2); % // compute squarederror           
        end

% // plot with optimal parameter
[T,Y] = ode45(@KIN1PARM, tspan, y0, [], k_opt);
figure
plot(ordinatus, abscissa,'ko', 'markersize',10,'markerfacecolor','black')
hold on
plot(T,Y, 'r--', 'linewidth', 2)


% // Another attempt with fminsearch and the integral form
t = ordinatus;
t_fit = linspace(min(ordinatus), max(ordinatus))
y = abscissa;

% create model function with parameters A0 = p(1) and k = p(2)
model = @(p, t) p(1)*exp(-p(2)*t);
e = @(p) sum((y - model(p, t)).^2); % minimize squared errors
p0 = [100, -0.1]; % an initial guess (positive A0 and probably negative k for exp. growth)
p_fit = fminsearch(e, p0); % Optimize 

% Add to plot
plot(t_fit, model(p_fit, t_fit), 'b-', 'linewidth', 2)

legend('location', 'best', 'data', 'ode45 with fixed y0', ...
    sprintf ('integral form: %5.1f*exp(-%.4f)', p_fit))
end

function dy = KIN1PARM(t,y,k)
%
% version : first order reaction
%   A  --> B
%   dA/dt = -k*A
%   integrated form A = A0*exp(-k*t)
%
dy = -k.*y; 
end 

结果如下所示。令我惊讶的是，y0 = 100 的初始猜测与找到的最佳A0 非常吻合。结果如下图：