Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I was writing a program to solve a system of equations with Runge-Kutta method, but it doesn't work well and I decided to check what Mathematica can say about it. It turns out that I have a stiff system of equations. The question is how can I show or prove that it is stiff?

Here is the code that I use:

kb = 10^-4;
Ch = 317;
gamma = 2350;
H = 0.2;
v = 50/3600;
C0 = 0;
z[R_] := Pi*R*H;
Wb[C_] := kb*Ch*(1 - C/Ch)^(4/3);
f[C_, R_] := (2*z[R]*Wb[C]*(1 - C/(2*Pi*gamma)) - (C - C0)*v)/(z[R]*R);
g[C_] := Wb[C]/(2*Pi*gamma);
s = First@
   NDSolve[{c'[t] == f[c[t], r[t]], r'[t] == g[c[t]], c[0] == 0, 
     r[0] == 0.15}, {c, r}, {t, 0, 25900000}, 
    StartingStepSize -> 0.01, Method -> "ExplicitRungeKutta"];
(*Plot[Evaluate[{c[t],r[t]}/.s],{t,0,25900000},PlotStyle->Automatic]*)
 {c[27.6], r[ 27.6], c[27.5], r[ 27.5]} /. s

And here is the error I got:

NDSolve::ndstf: At t == 27.532652604579695`, system appears to be stiff. Methods Automatic, BDF, or StiffnessSwitching may be more appropriate.

UPD: ODE's in TeX

$\Bigg\{ \begin{array}{} \dot{C}=\frac {2 (1 - \frac{C}{2 \gamma \pi}) W_b z - (C-C_0) \nu} {z(R) R} \\ \dot{R}=\frac{W_b}{2 \gamma \pi} \end{array}$

Where:

$W_b = k_bC_h(1-\frac{C}{C_h})^\frac{4}{3}$
$z=\pi R H$

and $\gamma$, $\nu$, $C_0$, $H$, $k_b$ and $C_h$ are constants.

share|improve this question
1  
There doesn't appear to be a precise definition: en.wikipedia.org/wiki/Stiff_equation –  Qiaochu Yuan May 16 '12 at 4:14
2  
It is qualitative, just like condition number for matrices. Do you have very different time scales in your ode? –  copper.hat May 16 '12 at 5:02
3  
In fact, one of the working definitions used for stiffness is "any (system of) equations where explicit methods like Runge-Kutta need to take very small step sizes (and thus become very inefficient) to maintain accuracy." –  J. M. May 16 '12 at 5:35
1  
Can you write out the equation in TeX using our built-in MathJax capabilities? I am not very good at parsing Mathematica in my head. –  Willie Wong May 16 '12 at 7:30
1  
Since your system is two-dimensional, you could visualize $(\dot C, \dot R)$ as a vector field over the $(C, R)$ space, which will give you a feel of the trajectories of your ODE. It looks like there is rapid descent in $C$ to a stable manifold near $\dot C = 0$, and then presumably much slower dynamics within that manifold; this is a common characteristic of stiff systems. (Also, you have a nasty singularity all along $R = 0$.) –  Rahul May 16 '12 at 20:36
show 2 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.