problem with time step in ODE ?!

3 次查看(过去 30 天)
21did21
21did21 2011-12-6
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:
  2 个评论
Walter Roberson
Walter Roberson 2011-12-6
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
21did21 2011-12-7
cste is a scalar, i have change the other error in my text.

请先登录,再进行评论。

回答(7 个)

Walter Roberson
Walter Roberson 2011-12-6
Have you tried a stiff solver such as ode23s ?

Matt Tearle
Matt Tearle 2011-12-6
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?

21did21
21did21 2011-12-6
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 ...
  1 个评论
Matt Tearle
Matt Tearle 2011-12-7
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.

请先登录,再进行评论。


21did21
21did21 2011-12-6
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...

Jan
Jan 2011-12-7
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.

21did21
21did21 2011-12-7
it's ok now i think, i have cleaned my code.
what do you think about it ?
  2 个评论
Walter Roberson
Walter Roberson 2011-12-9
I think that
Dyde theta + = const * (1 / y)-beta * y;
is still not valid code.
21did21
21did21 2011-12-10
yes it's true, the real thing it's this:
dYdE=theta+cste*(1/y)-beta*y;

请先登录,再进行评论。


21did21
21did21 2011-12-9
up :)

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by