ode45解決dif.equation與進一步適合exp.results

Question

我正在構建一個代碼來解決diff。方程：

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計算的值不同（對於該se應該大約爲7.0e-4 exp。數據）。

看着頌歌求解器的結果（SSE_minimization_1parm => val_theo）我發現ode函數給了我一個零向量。

有人可以幫助我，找出頌歌解算器的進展情況嗎？

非常感謝！

---

Answer 1

所以，現在我可以得到最好的。按照我的方式，我將座標值作爲時間，橫座標值表示爲您嘗試建模的測量值。此外，你似乎已經爲解算器設置了很多選項，我都省略了。首先使用ode45()使用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非常吻合。結果如下：