problem with time step in ODE ?!

Hi all,

I have an ODE to be solved but I have a concern with the use of this function, the following error message:

Warning: Failure at t = 5.844315e-003. Unable to Meet tolerances Without Reducing the integration step size below the Smallest value allowed (1.387779e-017) at time t.
> In ode23 at 362

I do not understand what is marked but I do not understand why or how to solve this ...

• This is my main script (I put that important):*

Global theta0;
Global beta0;
global const;
K1 = 0.5, K2 = 0.4;
theta = (a * m * K1) / 2;
beta = (K2 * a) / 2;
const = (b * M) / (2 * l);
[X Y] = ode23 (@ F10 [tempsI tempsF], initial);

Here is my function:

function Dyde=F10(t,y)
global theta;
global beta;
global const;
dYdE=theta+cste*(1/y)-beta*y;

for info here are the magnitude of my variables:

   initial=1e-05 it is arbitrarily because in theory it should be very close to 0
  tempsI= 0.0058
  tempsF = 0.0766
  const = 2.81e 18;
  theta = 3879;
  beta = 0.55;

I hope you can help me understand and solve my problem because I do not know what to do at all ...

thank you in advance: ccool:

[Walter Roberson]

The "= const" line in your main script seems to be missing a variable on the left hand side.

Is "const" a function or an array?

Your line "Dyde theta" in your F10 function is not valid syntax.

[21did21]

cste is a scalar, i have change the other error in my text.

Answer 1

[Walter Roberson]

Have you tried a stiff solver such as ode23s ?

Answer 2

[Matt Tearle]

I don't know exactly what you're trying to do, because your code is invalid in the last line of F10 and in your call to ode23. But if I do this:

const = 2.81e18;
theta = 3879;
beta = 0.55;
[X Y] = ode23s(@(t,y) F10(t,y,theta,beta,const), [0.0058,0.0766], 1e-5);

with

function Dyde = F10 (~, y, theta, beta, const)
Dyde = theta + const * (1 / y)-beta * y;

(which, by the way, is a much nicer way to handle extra parameters than using global), then I get the same error you get, mainly because you have a const/y term in your equation and const ~ 1e18 and y ~ 1e-5. So your ODE is just hideously scaled. Is there any way you can rescale/choose different units?

Answer 3

[21did21]

thanks a lot, for your answer !

=> Walter: i have tried ODE23s but it don't change something :(

=> Matt: i have change this 2.81e18 by 2.81e3 end now it's works!

i don't understand with with 2.81e18 cause problems...? there is no method to solve it ? if i change units const will be close to 1 but theta close to 1e18 ...

[Matt Tearle]

y is very small, but dy/de (or dy/dt or whatever it is) is huge. This means that even a very small timestep can lead to a big change in y, introducing all sorts of errors. This is just the nature of solving badly-scaled ODEs numerically.

What ode23/45/113 try to do is find a stepsize such that the resulting error is smaller than a given tolerance. For y ~ 1e-5, this default tolerance will be max(1e-8,1e-6) = 1e-6. To get an error that small when the derivative is ~ 1e23 will require a timestep small than the solver allows (~1e-17).

Perhaps you could try doing a nondimensionalization and see what parameters drop out and what sizes they are likely to be. If there's no way around having a wide range of scales, you might need to think about a different approach.

Answer 4

[21did21]

perhaps i can add some options but i don't know how i can do this because in the help we don't have example for ODE...

Answer 5

[Jan]

This is no proper code:

Global theta0;
Global beta0;
global const;
K1 = 0.5, K2 = 0.4;
theta = (a * m * K1) / 2;
beta = (K2 * a) / 2;
= const (b * M) / (2 * l);
[X {1} {1} Y] = ode23 (@ F10 [tempsI tempsF {1} {1}], initial {1});

It is not possible to guess, what's going wrong, if we start from this piece on non-Matlab code. Please fix this by editing the question.

Answer 6

[21did21]

it's ok now i think, i have cleaned my code.

what do you think about it ?

[Walter Roberson]

I think that

Dyde theta + = const * (1 / y)-beta * y;

is still not valid code.

[21did21]

yes it's true, the real thing it's this:

dYdE=theta+cste*(1/y)-beta*y;

Answer 7

[21did21]

up :)