Parameter estimation using fminsearch

Hi,

I'm trying to determine 5 unknown parameters from a ODE system. But I get this warning which I don't understand.

-----------------------------------------------------------------

Warning: Failure at t=7.437938e+000. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.421085e-014) at time t. > In ode15s at 819 In trial_least at 41 In fminsearch at 326 In Main at 28 ??? Error using ==> deval at 143 Attempting to evaluate the solution outside the interval [7.000000e+000, 7.437938e+000] where it is defined.

Error in ==> trial_least at 51 Data_cal_1 = deval(Soln, Time_1);

Error in ==> fminsearch at 326 x(:) = xe; fxe = funfcn(x,varargin{:});

Error in ==> Main at 28 [ParamFinal, Func] = fminsearch(@trial_least, Param, options);

----------------------------------------------------------------

Here's my coding:

------------------------------------------------------------- Main.m -------------------------------------------------------------

clear all

%%% the initial guess for fminsearch

beta = 0.1;

gamma = 0.04;

k = 3;

L = -20;

D_0 = 0.1;

Param = [beta; gamma; k; L; D_0];

Func = Least_Squares(Param);

options = optimset('Display','Final', 'MaxIter', 10000, 'MaxFunEvals', ... 10000, 'TolX',1e-8,'TolFun',1e-8);

[ParamFinal, Func] = fminsearch(@Least_Squares, Param, options);

fprintf('%s \n' , 'Final Optimum Parameters:')

ParamFinal

----------------------------------------------------------------- Least_Squares.m ------------------------------------------------------------------

function [SSE_val] = Least_Squares(Param)

global Data_1 Data_2 Data_3 Data_4 Data_5;

%% Actual data

Data_1 = [0.125 0.1875 0.25 0.25 0.25];

Data_2 = [0.0625 0.09375 0.125 0.125 0.125];

Data_3 = [0.125 0.1875 0.25 0.3125 0.375];

Data_4 = [0.3125 0.40625 0.5 0.5 0.5];

Data_5 = [0.5 0.5625 0.625 0.625 0.625];

%% define the parameter vector

beta = Param(1);

gamma = Param(2);

k = Param(3);

L = Param(4);

D_0 = Param(5);

%% the integration step

dt = 0.01;

%% t_p

t_p_1 = 7;

t_p_2 = 14;

t_p_3 = 21;

t_p_4 = 35;

t_p_5 = 0;

%% the initial time

t_i_1 = t_p_1;

t_i_2 = t_p_2;

t_i_3 = t_p_3;

t_i_4 = t_p_4;

t_i_5 = t_p_5;

%% the final time

t_f = 150;

%% the integration time

tspan_1 = t_i_1 : dt : t_f;

tspan_2 = t_i_2 : dt : t_f;

tspan_3 = t_i_3 : dt : t_f;

tspan_4 = t_i_4 : dt : t_f;

tspan_5 = t_i_5 : dt : t_f;

%% the initial time

y0 = [0 0];

options = odeset ('RelTol', 1e-6, 'AbsTol', [1e-6 1e-6]);

IndRes_1 = @(t,y)IR_1 (t, y, beta, gamma, k, L, D_0);

IndRes_2 = @(t,y)IR_2 (t, y, beta, gamma, D_0);

Soln_1 = ode15s(IndRes_1, tspan_1, y0, options);

Soln_2 = ode15s(IndRes_1, tspan_2, y0, options);

Soln_3 = ode15s(IndRes_1, tspan_3, y0, options);

Soln_4 = ode15s(IndRes_1, tspan_4, y0, options);

Soln_5 = ode15s(IndRes_2, tspan_5, y0, options);

%% time to record the disease

%% for t_p = 7

t1_d1 = t_p_1 + 12;

t1_d2 = t_p_1 + 16;

t1_d3 = t_p_1 + 20;

t1_d4 = t_p_1 + 24;

t1_d5 = t_p_1 + 27;

%% for t_p = 14

t2_d1 = t_p_2 + 12;

t2_d2 = t_p_2 + 16;

t2_d3 = t_p_2 + 20;

t2_d4 = t_p_2 + 24;

t2_d5 = t_p_2 + 27;

%% for t_p = 21

t3_d1 = t_p_3 + 12;

t3_d2 = t_p_3 + 16;

t3_d3 = t_p_3 + 20;

t3_d4 = t_p_3 + 24;

t3_d5 = t_p_3 + 27;

%% for t_p = 35

t4_d1 = t_p_4 + 12;

t4_d2 = t_p_4 + 16;

t4_d3 = t_p_4 + 20;

t4_d4 = t_p_4 + 24;

t4_d5 = t_p_4 + 27;

%% for t_p = 0

t5_d1 = t_p_5 + 12;

t5_d2 = t_p_5 + 16;

t5_d3 = t_p_5 + 20;

t5_d4 = t_p_5 + 24;

t5_d5 = t_p_5 + 27;

%% calling solution at the specific time %% for t_p = 7

Time_1 = [t1_d1 t1_d2 t1_d3 t1_d4 t1_d5];

Data_cal_1 = deval(Soln_1, Time_1);

Data_D_1 = Data_cal_1(2, :);

%% for t_p = 14

Time_2 = [t2_d1 t2_d2 t2_d3 t2_d4 t2_d5];

Data_cal_2 = deval(Soln_2, Time_2);

Data_D_2 = Data_cal_2(2, :);

%% for t_p = 21

Time_3 = [t3_d1 t3_d2 t3_d3 t3_d4 t3_d5];

Data_cal_3 = deval(Soln_3, Time_3);

Data_D_3 = Data_cal_3(2, :);

%% for t_p = 35

Time_4 = [t4_d1 t4_d2 t4_d3 t4_d4 t4_d5];

Data_cal_4 = deval(Soln_4, Time_4);

Data_D_4 = Data_cal_4(2, :);

%% for t_p = 0

Time_5 = [t5_d1 t5_d2 t5_d3 t5_d4 t5_d5];

Data_cal_5 = deval(Soln_5, Time_5);

Data_D_5 = Data_cal_5(2, :);

%% calculate SSE %% this is the calculate data

Data_D = [Data_D_1; Data_D_2; Data_D_3; Data_D_4; Data_D_5];

%% this is the actual data

Data = [Data_1; Data_2; Data_3; Data_4; Data_5];

SSE = 0;

for i = 1 : 5

    for x = 1 : 5

        SSE = SSE + ((Data_D(i,x) + D_0) - Data(i,x))^2;

    end

end

SSE_val = SSE;

------------------------------------------------------------------ IR_1.m -----------------------------------------------------

function dy = IR_1 (t, y, beta, gamma, k, L, D_0)

dy = [(k*t/ (t^2 + L))* (1 - y(1) - y(2)- D_0) - gamma * y(1); ... beta * (y(2) + D_0)* (1 - y(1) - y(2) - D_0)];

--------------------------------------------------------------- IR_2.m --------------------------------------------------------------- function dy = IR_2 (t, y, beta, gamma , D_0)

dy = [-gamma * y(1); beta * (y(2) + D_0)* (1 - y(1) - y(2) - D_0)];

Really appreciate your suggestions / comments. Thanks!

Syaza