Re: How to solve a stiff ODE system?

Orlando wrote:
Narasimham wrote:
Orlando wrote:
d Y
---- = f(t, Y), Y(t_0) = Y_0
d t
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.

Suspect stiffness is not a problem,cannot say further as you do not
indicate f(t,Y).For eps =10^-8, if initial condition Y(t_eps) = Y_0 is
used and ode behaves well, then problem was with usage of first order
equation.To fix such cases switch over to a second order differential
equation providing two initial conditions.

Thank you very much for your answer!

The ode system of mine is as follows:

dy(1) / dt = y(2)
dy(2) / dt = - (y(3) - y(4)) / epsilon0 / E^2
dy(3) / dt = 1 / epsilon0 / y(2) / E^2 * (y(3)^2 - y(3) * y(4) * (1 -
epsilon0 * Ri / mu_+ / e))
dy(4) / dt = 1 / epsilon0 / y(2) / E^2 * (y(3) * y(4) * (1 - epsilon0 *
Ri / mu_- / e) - y(4)^2)

Where,
epsilon0: permittivity of vacuum 8.85e-12
E = E(t): the electric field. it is a function of t. After integrating
one step of the ODE, the electric field should be recalculated from
y(1), y(2), y(3), y(4), and this formula doesn't belong to the ODE
system, so it is not given here.
Ri: the recombination coefficient of bipolar charge. 2.2e-12
mu_+, mu_-: the mobility of positive and negative charge, 1.6e-4
e: elementary charge 1.6e-19

The above system is used to model the space charge flow of the high
voltage direct current power line.

I still don't quite follow you: how should I switch to a second order
differential equation? What is the reason of doing this? And how to
apply this procedure to my ODE system?

Thank you again!!

Before going to second order, at first redefine terms so that you
compute quantities in absolute value range (.01-10).Tried to solve the
ODE numerically.It crashes at start of any chosen time interval, even
as an initial value problem,found no errors. Perhaps you need to at
first check the physics modeling even as an IVP.In a simulation program
start with the simplest electrostatic situations of equipotentials or
force lines and equations depicting them.When it works,gradually
increase its complexity in steps.That traps the error at a particular
stage.

Narasimham

.

Relevant Pages

  • Re: How to solve a stiff ODE system?
    ... Unable to meet integration ... The ode system of mine is as follows: ... the recombination coefficient of bipolar charge. ... apply this procedure to my ODE system? ...
    (sci.math)
  • Re: How to solve a stiff ODE system?
    ... Unable to meet integration ... the recombination coefficient of bipolar charge. ... I still don't quite follow you: how should I switch to a second order ... apply this procedure to my ODE system? ...
    (sci.math)
  • Re: How to deal with such ODE with fortran?
    ... I have to an integrate different ODE for each kp. ... without causing confusion in the recursive procedure. ... RHscalls P, but the subroutine doing integration only ... I am not familiar with this kind of equation, so I am now writing codes ...
    (comp.lang.fortran)
  • Re: Runge-Kutta question
    ... point of phase transitions) and, therefore, a numerical solution goes crazy. ... Try regularizing the heat capacity by replacing the jump with a steep, ... Your ODE will then have an analytical ... the integration step to really small values. ...
    (sci.math.num-analysis)
  • Re: How to deal with such ODE with fortran?
    ... I have to an integrate different ODE for each kp. ... without causing confusion in the recursive procedure. ... RHscalls P, but the subroutine doing integration only ... available literature for suitable numerical methods? ...
    (comp.lang.fortran)