Problems with ode23s

I have this warning on matlab

Warning: Failure at t=4.447832e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.580187e-14) at time t.

the important code is this and i dont know how to solve it:

t0 = 0.001;
tf = 144;
Opciones = odeset('RelTol', 1e-4, 'AbsTol', 1e-7);
[t,X] = ode23s('SedimenEstNN', [t0,tf], [X(1) X(2) X(3) X(4)...
    X(5) X(6) X(7) X(8) X(9) X(10) X(11) X(12) X(13) X(14)...
    X(15) X(16) X(17)], Opciones);

function dX=SedimenEstNN(t,X)

global A H K N P Va be 

U = SedimenEntNN(t);
Qd = P*U(1);
Eda = Va*U(2)/(24*Qd*X(16));
dX(17) = Eda;
Edad =X(17)/t;
%more code...

function U=SedimenEntNN(t)

global A H K N P Va be 

A = 200;     H = 2;     K = 7e-13;
N = 6;       be = 2.2;  Va = 600;

w = 2*pi/24*(t-10);
R = 0.75;
P = 0.0156;

if t>118
    U(1) = 100;   U(2) = 2500;   U(3) = R*U(1);
else
    U(1) = 100 + 50*sin(w);         
    U(2) = 2500 + 500*sin(w);       
    U(3) = R*U(1);                  
end

thank you before all and if i need to put more code tell me please

[Jan]

 [X(1) X(2) X(3) X(4)...
    X(5) X(6) X(7) X(8) X(9) X(10) X(11) X(12) X(13) X(14)...
    X(15) X(16) X(17)]

?? What about:

X(1:17)

Answer 1

[Torsten]

It's hard to tell why ode23s fails for your problem.

Some options how to proceed:

1. Make sure that you transposed the dX-array at the end of function "SedimenEstNN".

2. The vector of initial guesses [X(1) X(2) X(3) X(4) X(5) X(6) X(7) X(8) X(9) X(10) X(11) X(12) X(13) X(14) X(15) X(16) X(17)] should be a column vector, not a row vector. Did you give values to the X(i)'s before you called ode23s ?

3. Check the solution up to the time where the solver fails.

4. Check the dX values at the time where the solver fails.

5. Choose smaller tolerances.

6. Try a stiff solver (e.g. ode15s).

Good luck.

Best wishes

Torsten.

[Pablo R]

Hi Torsten,

i tried all of that but still having the same problem. i show you by the way,

1)

dX = [dX(1); dX(2); dX(3); dX(4); dX(5); dX(6); dX(7); dX(8);...
    dX(9); dX(10); dX(11); dX(12); dX(13); dX(14); dX(15);...
    dX(16); dX(17)];
2)

X(1) = 100;      X(2) = 75;       X(3) = 2;        X(4) = 2500;
X(5) = 100;      X(6) = 75;       X(7) = 3;        X(8) = 100;
X(9) = 75;       X(10) = 3;       X(11) = 2500;    X(12) = 2500;
X(13) = 223;     X(14) = 175;     X(15) = 30;      X(16) = 6000;
X(17) = 1e-3;

5) same problem

6) same problem

i will see the values at that time

but thank you so much for your answer!

Answer 2

[Jan]

 [X(1) X(2) X(3) X(4)...
    X(5) X(6) X(7) X(8) X(9) X(10) X(11) X(12) X(13) X(14)...
    X(15) X(16) X(17)]

?? What about:

X(1:17)

Do not provide the function to be integrated as string 'SedimenEstNN', but as function handle @SedimenEstNN. This is the preferred method since Matlab 6.5 (over 15 years now...).

The function to be integrated is not smooth at t=118. This is a really bad idea and Matlab's integrators are not designed to handle such jumps. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047

Use the reached time as endpoint and look at the trajectory:

t0 = 0.001;
tf = 4.447832e+00 - 10 * eps(4.447832);
Opciones = odeset('RelTol', 1e-4, 'AbsTol', 1e-7);
[t,X] = ode23s('SedimenEstNN', [t0,tf], [X(1) X(2) X(3) X(4)...
  X(5) X(6) X(7) X(8) X(9) X(10) X(11) X(12) X(13) X(14)...
  X(15) X(16) X(17)], Opciones);

plot(t, X);   % Or X.'? I cannot remember this detail

Do you see a pole? If so, it is a mathematical problem and not a fault of the integrator.