A basic ODE where Runge Kutta doesn't work?

  • Context: Graduate 
  • Thread starter Thread starter h1a8
  • Start date Start date
  • Tags Tags
    Ode Runge kutta Work
Join the discussion
Registration is free. Ask a follow-up in this thread, or start your own.
6 replies · 5K views
h1a8
Messages
86
Reaction score
4
From a cubic function where y(0)=1, y(1)=0, and where there is a local max at y(5/13) I created a basic separable differential equation problem. I wanted to analyze how well different ordered Runge Kutta methods works in an interval [0,1]. Here it is:

dy/dt=-6(6/13)1/3(y-343/468)2/3 , y(0)=1

This ODE yields the cubic solution of

y=1/468(-12t+5)3+343/468

Now it is clear that y(1)=0

But using the several Runge Kutta programs with various computer software (mathematica, ti-nspire cas, mathstudio, etc.) yields a complex solutions for y(1). For example, using the classical RK4 with h=.1 yields
y(1)=0.718779+.005811i.

I don't see how the programs get a complex solution when all the functions have no even roots. Does anyone what is going on with these programs?
 
Physics news on Phys.org
You have negative numbers in that root for y<343/468. All of the mentioned programs are won't to go to complex numbers ASAP. Those programs will take the principal root. For example, the principal root of (-1)2/3 is (-1+i√3)/2, not 1.
 
Last edited:
Yep.

Try dy/dt=-6 (6/13)1/3 ((y-343/468)(y-343/468))1/3 instead.

By multiplying (y-343/468) with itself first, you make sure that the argument to the power (1/3) is non-negative.
This way it will do what you intended.
 
I like Serena said:
Yep.

Try dy/dt=-6 (6/13)1/3 ((y-343/468)(y-343/468))1/3 instead.

By multiplying (y-343/468) with itself first, you make sure that the argument to the power (1/3) is non-negative.
This way it will do what you intended.

Thanks that works. I tried ((y-343/468)2)1/3 and that works too.
Didn't know the programs calculated the principle cube root of negatives by default even if they were complex.
 
Well, basically they don't have a choice.

What do you think (-1)^0.667 is?
 
I like Serena said:
Try dy/dt=-6 (6/13)1/3 ((y-343/468)(y-343/468))1/3 instead.
Nice trick.

h1a8, this works because of the identity (ab)c=abc.

This identity is valid if a is positive and b and c are real. It is not valid in general. (The same goes with a lot of other exponentiation identities.)
 
Last edited: