| |
 |
|
|
Science Forum Index » Math - Numerical Analysis Forum » How to solve a stiff ODE system?
Page 1 of 1
|
| Author |
Message |
| Orlando |
 Posted: Fri Dec 01, 2006 11:47 pm |
|
|
|
Guest
|
Hello everybody!
I want to solve a system of first order ODEs with a given initial
value:
d Y
---- = f(t, Y) Y(t_0) = Y_0
d t
Where, f(t, Y) is nonlinear.
Up till now, I've tried to deal with this IVP using a series of ode**
solvers (stiff or nonstiff) in Matlab. However, when integrating to
some certain time points, a warning message is thrown out from the
ode** solver as follows:
Warning: Failure at t=1.533227e-003. Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (5.447117e-018) at time t.
So I examined the Jacobian matrix (\frac{\partial f}{\partial Y}) of
the nonlinear mapping at the current time point and obtained its
eigenvalues as follows:
lambda1 = -1.72864992567911e9
lambda2 = 0.29658911109064e9
lambda3 = -0.17145123351451e9
lambda4 = 0
We can see that this Jacobian matrix has a negative eigenvalue with
very large absolute value, and also a very large positive eigenvalue.
According to the definition of the stiffness of the ODEs, which is
given by Shampine and Gear: a ODE system is stiff, when:
1) the real parts of all the eigenvalues of the Jacobian matrix should
not be very large positive real numbers;
2) there should be at least one eigenvalue of the Jacobian matrix,
which has a negative real part with very large absolute value;
3) ...
My present situation complies with second requirement of the stiffness
definition, but not the first one, due to the existence of the
eigenvalue: 0.29658911109064e9.
Then my questions are:
1) Is my ode system stiff?
2) if so, why do the stiff ode solvers of Matlab failed at my problem?
What will be the possible reasons?
3) if my ode system is not stiff, then is the method for solving stiff
ode also suitable for my problem? And if the method of solving stiff
ode does not work, are there any other methods that can be effective?
I am very grateful for any suggestions and opinions, thanks a lot for
your help! |
|
|
| Back to top |
|
| Fred Krogh |
| Posted: Sun Dec 03, 2006 1:59 pm |
|
|
|
Guest
|
Orlando wrote:
Quote: Hello everybody!
I want to solve a system of first order ODEs with a given initial
value:
d Y
---- = f(t, Y) Y(t_0) = Y_0
d t
Where, f(t, Y) is nonlinear.
Up till now, I've tried to deal with this IVP using a series of ode**
solvers (stiff or nonstiff) in Matlab. However, when integrating to
some certain time points, a warning message is thrown out from the
ode** solver as follows:
Warning: Failure at t=1.533227e-003. Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (5.447117e-018) at time t.
So I examined the Jacobian matrix (\frac{\partial f}{\partial Y}) of
the nonlinear mapping at the current time point and obtained its
eigenvalues as follows:
lambda1 = -1.72864992567911e9
lambda2 = 0.29658911109064e9
lambda3 = -0.17145123351451e9
lambda4 = 0
We can see that this Jacobian matrix has a negative eigenvalue with
very large absolute value, and also a very large positive eigenvalue.
According to the definition of the stiffness of the ODEs, which is
given by Shampine and Gear: a ODE system is stiff, when:
1) the real parts of all the eigenvalues of the Jacobian matrix should
not be very large positive real numbers;
2) there should be at least one eigenvalue of the Jacobian matrix,
which has a negative real part with very large absolute value;
3) ...
My present situation complies with second requirement of the stiffness
definition, but not the first one, due to the existence of the
eigenvalue: 0.29658911109064e9.
Then my questions are:
1) Is my ode system stiff?
2) if so, why do the stiff ode solvers of Matlab failed at my problem?
What will be the possible reasons?
3) if my ode system is not stiff, then is the method for solving stiff
ode also suitable for my problem? And if the method of solving stiff
ode does not work, are there any other methods that can be effective?
I am very grateful for any suggestions and opinions, thanks a lot for
your help!
There are a number of things to consider. If you are interested in the
part of the solution brought in by the positive eigenvalue and are using
a relative error test then your system is at best only slightly stiff.
(The large positive eigenvalue limits the stepsize very nearly as much
at the negative ones.) If you are using an absolute error test then the
problem you are having is likely due to the positive eigenvalue causing
the solution to blow up, and not to stiffness. (I'm guessing that this
is the source of your problem.) Is there a chance that you should be
solving this problem as a boundary value problem?
I favor a definition of stiffness which depends on the type of method
that works best. With this definition, no initial value problem is
stiff when integrated on a sufficiently short interval or if a positive
eigenvalue has magnitude greater than the most negative one.
Finally, it somewhat interesting that the BDF based methods will solve
your type of problem if you want to damp out the part of the solution
due to the positive eigenvalue if you always use a sufficiently LARGE
stepsize.
Fred |
|
|
| Back to top |
|
| Peter Spellucci |
| Posted: Mon Dec 04, 2006 1:54 am |
|
|
|
Guest
|
In article <1165031273.028886.50790@n67g2000cwd.googlegroups.com>,
"Orlando" <tianjh04@gmail.com> writes:
Quote: Hello everybody!
I want to solve a system of first order ODEs with a given initial
value:
d Y
---- = f(t, Y) Y(t_0) = Y_0
d t
Where, f(t, Y) is nonlinear.
Up till now, I've tried to deal with this IVP using a series of ode**
solvers (stiff or nonstiff) in Matlab. However, when integrating to
some certain time points, a warning message is thrown out from the
ode** solver as follows:
Warning: Failure at t=1.533227e-003. Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (5.447117e-018) at time t.
So I examined the Jacobian matrix (\frac{\partial f}{\partial Y}) of
the nonlinear mapping at the current time point and obtained its
eigenvalues as follows:
lambda1 = -1.72864992567911e9
lambda2 = 0.29658911109064e9
lambda3 = -0.17145123351451e9
lambda4 = 0
We can see that this Jacobian matrix has a negative eigenvalue with
very large absolute value, and also a very large positive eigenvalue.
According to the definition of the stiffness of the ODEs, which is
given by Shampine and Gear: a ODE system is stiff, when:
1) the real parts of all the eigenvalues of the Jacobian matrix should
not be very large positive real numbers;
2) there should be at least one eigenvalue of the Jacobian matrix,
which has a negative real part with very large absolute value;
3) ...
My present situation complies with second requirement of the stiffness
definition, but not the first one, due to the existence of the
eigenvalue: 0.29658911109064e9.
Then my questions are:
1) Is my ode system stiff?
2) if so, why do the stiff ode solvers of Matlab failed at my problem?
What will be the possible reasons?
3) if my ode system is not stiff, then is the method for solving stiff
ode also suitable for my problem? And if the method of solving stiff
ode does not work, are there any other methods that can be effective?
I am very grateful for any suggestions and opinions, thanks a lot for
your help!
this would not be a stiff equation: your solution would be totally dominated by
the extremely fast growing part, and all solutions in the neighborhood of this
specific one would also grow so rapidly. I guess you are running into a
singularity . using merely a reltol criterion instead of an abstol would also not
help, because the next thing would be an overflow and you are computing inf's only
you should examine the problem (or your code?)
hth
peter |
|
|
| Back to top |
|
| bv |
| Posted: Thu Dec 07, 2006 10:35 pm |
|
|
|
Guest
|
Peter Spellucci wrote:
Quote:
So I examined the Jacobian matrix (\frac{\partial f}{\partial Y}) of
the nonlinear mapping at the current time point and obtained its
eigenvalues as follows:
lambda1 = -1.72864992567911e9
lambda2 = 0.29658911109064e9
lambda3 = -0.17145123351451e9
lambda4 = 0
this would not be a stiff equation: your solution would be totally dominated by
the extremely fast growing part, and all solutions in the neighborhood of this
specific one would also grow so rapidly...
Professor, you're way off today, so much so that your answer should earn
you an "F" or a loss of tenure. In fact, it's as bad as some of the
opinions you gave in the "Merson integration" thread - end of sermon.
Now, you would be right had it not been for the (slow) mode associated
with the 4th eigenvalue. Blindly applying a stiffness ratio definition
never was a good idea - I'm speculating that's what you did and dam the
numerical evidence to the contrary. For your convenience, a couple of
definitions that can be easily understood.
wikipedia: "a stiff equation is a differential equation for which
certain numerical methods for solving the equation are numerically
unstable, unless the step size is taken to be extremely small..."
elsewhere: "stiff differential equations are those with two or more
widely differing scales. For example, the solution could have a
component that quickly becomes insignificant, and another component that
changes much more slowly."
--
sdx
http://www.sdynamix.com |
|
|
| Back to top |
|
| DP |
| Posted: Mon Dec 11, 2006 4:14 am |
|
|
|
Guest
|
Orlando wrote:
Quote: Hello everybody!
I want to solve a system of first order ODEs with a given initial
value:
d Y
---- = f(t, Y) Y(t_0) = Y_0
d t
Where, f(t, Y) is nonlinear.
Up till now, I've tried to deal with this IVP using a series of ode**
solvers (stiff or nonstiff) in Matlab. However, when integrating to
some certain time points, a warning message is thrown out from the
ode** solver as follows:
Warning: Failure at t=1.533227e-003. Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (5.447117e-018) at time t.
So I examined the Jacobian matrix (\frac{\partial f}{\partial Y}) of
the nonlinear mapping at the current time point and obtained its
eigenvalues as follows:
lambda1 = -1.72864992567911e9
lambda2 = 0.29658911109064e9
lambda3 = -0.17145123351451e9
lambda4 = 0
We can see that this Jacobian matrix has a negative eigenvalue with
very large absolute value, and also a very large positive eigenvalue.
According to the definition of the stiffness of the ODEs, which is
given by Shampine and Gear: a ODE system is stiff, when:
1) the real parts of all the eigenvalues of the Jacobian matrix should
not be very large positive real numbers;
2) there should be at least one eigenvalue of the Jacobian matrix,
which has a negative real part with very large absolute value;
3) ...
My present situation complies with second requirement of the stiffness
definition, but not the first one, due to the existence of the
eigenvalue: 0.29658911109064e9.
Then my questions are:
1) Is my ode system stiff?
2) if so, why do the stiff ode solvers of Matlab failed at my problem?
What will be the possible reasons?
3) if my ode system is not stiff, then is the method for solving stiff
ode also suitable for my problem? And if the method of solving stiff
ode does not work, are there any other methods that can be effective?
I am very grateful for any suggestions and opinions, thanks a lot for
your help!
Did you try to change of variables? If you would show your system
readers could provide ideas along this line. Sometimes singularities
can be removed, or stiff systems transformed into regular ones.
DP |
|
|
| Back to top |
|
| |
|
Page 1 of 1
All times are GMT - 5 Hours
The time now is Mon Dec 01, 2008 10:51 am
|
|