Message Boards Message Boards

0
|
9984 Views
|
6 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Runge-Kutta Method, stiffness occur, how to solve it?

Posted 10 years ago

Hello everyone, I'm a new user to Mathematica and i have been trying to solve my equations using Explicit Runge-Kutta Method. The problem is that I have encountered "stiffness" problem, which shouldn't happen according to the calculation.

At first, I need to solve this equation by using 4th Order Runge-Kutta Method, but, it seems like Mathematica have a build it package (Explicit Runge-Kutta Method). I tried to use it and problem occur. Can anybody advice on what is the problem?, Should it be my own equation missed out something?

POSTED BY: Thai Kee Gan
6 Replies
Posted 10 years ago

Hello guys,

I played around with NDSolve and I found out that ExplicitEuler Method works instead of ExplicitRungeKutta Method Does ExplicitEuler Method not encounter "infinity" as stiffness?

Doesn't Work NDSolve[{y'[x] == y[x]^3 (1 - y[x]) x^2 - (1 - y[x])^3 x^2 + 3/8 y[x] (1/2 y[x]^(-1/2)), y[0] == 0.}, y, {x, 0, 1}, Method -> "ExplicitRungeKutta", "StartingStepSize" -> 1/1000]

Work NDSolve[{y'[x] == y[x]^3 (1 - y[x]) x^2 - (1 - y[x])^3 x^2 + 3/8 y[x] (1/2 y[x]^(-1/2)), y[0] == 0.}, y, {x, 0, 1}, Method -> "ExplicitEuler", "StartingStepSize" -> 1/10]

NDSolve[{y'[x] == y[x]^3 (1 - y[x]) x^2 - (1 - y[x])^3 x^2 + 3/8 y[x] (1/2 y[x]^(-1/2)), y[0] == 0.}, y, {x, 0, 1}, Method -> "ExplicitEuler", "StartingStepSize" -> 1/1000]

Btw, I get 2 different graphs for different step size.

Also, I would like to know how to input my equation into the Runge Kutta Plug-in Method. It was in the tutorial guide under ExplicitRungeKuttaMethod but what was written is just the basic formula.

ClassicalRungeKutta[__]["Step"[f, t, h, y, yp]] := Block[{deltay, k1, k2, k3, k4}, k1 = yp; k2 = f[t + 1/2 h, y + 1/2 h k1]; k3 = f[t + 1/2 h, y + 1/2 h k2]; k4 = f[t + h, y + h k3]; deltay = h (1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4); {h, deltay} ];

What does this command means? ["Step"[f, t, h, y, yp_]]

Thanks your the help :D

POSTED BY: Thai Kee Gan
Posted 10 years ago

Thank you Kay, I'd assume that y[x]^(-1/2) will be automatically sorted out when I run command. I guess that was my mistake. Also, how do I combine NDSolve and Piecewise together?

For example, y'[x]==Piecewise[{{x^2 (-1 + y[x])^3 + (3 Sqrt[y[x]])/16 - x^2 (-1 + y[x]) y[x]^3, 0 <= x <=0.5}, {x^2 (-1 + y[x])^3 + (3 Sqrt[y[x]])/16 - x^2 (-1 + y[x]) y[x]^3+30, 0.5 <= x <=1}}]

I tried NDSolve[{y'[x] == Piecewise[{{x^2 (-1 + y[x])^3 + (3 Sqrt[y[x]])/16 - x^2 (-1 + y[x]) y[x]^3, 0 <= x <= 0.5}, {x^2 (-1 + y[x])^3 + (3 Sqrt[y[x]])/16 - x^2 (-1 + y[x]) y[x]^3 + 30, 0.5 <= x <= 1}}], y[0] == 0}, y, {x, 0.000, 1}, Method -> Automatic]

But I get nothing when I plot the graph.

POSTED BY: Thai Kee Gan

Your function has a square root term of negative value:

In[182]:= Clear["Global`*"]

In[183]:= Simplify[
 y[x]^3 (1 - y[x]) x^2 - (1 - y[x])^3 x^2 + 3/8 y[x] (1/2 y[x]^(-1/2))]

Out[183]= 
x^2 (-1 + y[x])^3 + (3 Sqrt[y[x]])/16 - x^2 (-1 + y[x]) y[x]^3

In[184]:= s = 
 NDSolve[{y'[x] == 
    x^2 (-1 + y[x])^3 + (3 Sqrt[Abs[y[x]]])/16 - 
     x^2 (-1 + y[x]) y[x]^3, y[0] == 0}, y, {x, 0.000, 1}, 
  Method -> Automatic]

Out[184]= {{y -> InterpolatingFunction[{{0., 1.}}, <>]}}

It works if I take absolute value underneath the square root. I also used Automatic as the integration method.

POSTED BY: Kay Herbert
Posted 10 years ago

Hello guys, thanks for the replies. My equation is:

NDSolve[{y'[x] == 
   y[x]^3 (1 - y[x]) x^2 - (1 - y[x])^3  x^2 + 
    3/8 y[x] (1/2 y[x]^(-1/2)), y[0] == 0.`}, y, {x, 0, 1`}, 
 Method -> "ExplicitRungeKutta", "StartingStepSize" -> 1/1000]

Is there any mistake on how I load the Runge-Kutta Package? Or is there any mistake in the equation?

POSTED BY: Thai Kee Gan

May I make explicit a point which was implicit already in David's answer: The 'normal' way to obtain numerical solutions for ODEs from Mathematica is to call NDSolve and leave it to the system to find out among the many built-in methods and merthod variants the one which seems to be most suited for your equation. Since it definitively has methods which are designed to cope with stiff equations it should come up with a reliable solution whenever it exists mathematically. Of course, it can only make you learn and understand if you compare the behavior of standard methods for cases in which some of them get into trouble. Unfortunately the documentation on the 'Method' option of NDSolve is poor, so that some experimentation is needed to see which methods are implemented.

POSTED BY: Ulrich Mutze

Could you include the system of equations you are trying to solve? In some recent work, I have had 'stiffness' warning messages from NDSolve for certain sets of input parameters when the same set of equations solve perfectly well under different initial conditions.

POSTED BY: David Mackay
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract