A basic ODE where Runge Kutta doesn't work?

[h1a8]

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)³+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?

 

[D H]

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.

 

[I like Serena]

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.

 

[h1a8]

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)²)^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.

 

[I like Serena]

Well, basically they don't have a choice.

What do you think (-1)^0.667 is?

 

[D H]

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 (a^b)^c=a^bc.

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.)

 

[I like Serena]

To illustrate the problems of fractional powers of negative numbers, here's a nice example for real numbers:

−27 = (−27)^{((2/3)⋅(3/2))} = ((−27)^2/3)^3/2 = 9^3/2 = 27


Taken from: http://en.wikipedia.org/wiki/Exponentiation#Rational_powers.