6
$\begingroup$

Bug noticed for Mathematica 8 and persisting through 11.3; Submitted to Wolfram, Inc. as CASE:3729158

I want to solve the following non-linear first order difference equation:

$$ k_{t+1} = (1-a) k_t^{a} $$

with $0 < a < 1$, and $k_0 = \text{const.}$ .

Of course, the solution is easily found, but I want Mathematica to solve it.

Using

RSolve[{ k[t + 1] == (1-a) k[t]^a, k[0] == k0}, k[t], t ]

gives

{{k[t] -> E^(((-1 + a^t) (I π + Log[-1 + a]))/(-1 + a)) k0^a^t}}

The solution is much simpler when considering the restriction on $a$. How can I tell Mathematica that $a$ is between $0$ and $1$ when applying the RSolve command?

P.S. I've played around with the Reduce command, using assumptions, but that did not help.

$\endgroup$
  • 1
    $\begingroup$ Would someone verify that there is a bug with solution returned by RSolve? $\endgroup$ – Michael E2 Dec 23 '13 at 1:14
  • 1
    $\begingroup$ @MichaelE2 Bug persists in 11.0.1. Do you know whether anyone has reported in to Wolfram, Inc.? $\endgroup$ – bbgodfrey Oct 1 '16 at 21:19
  • 1
    $\begingroup$ @bbgodfrey No, I don't remember or know. $\endgroup$ – Michael E2 Oct 2 '16 at 2:56
  • 1
    $\begingroup$ @MichaelE2 I have submitted it. $\endgroup$ – bbgodfrey Oct 2 '16 at 3:46
4
$\begingroup$

Well, the answer appears to be wrong:

rsol = First @ RSolve[{k[t + 1] == (1 - a) k[t]^a, k[0] == k0}, k[t], t];
Table[k[t] /. rsol /. {a -> 1/2, k0 -> 1}, {t, 0, 5}] // N // Chop

(* {1., 0.5, -0.353553, 0. - 0.297302 I, 0.192776 - 0.192776 I, 0.241196 - 0.0999066 I} *)

All the terms in the recurrence relation are positive, but we're get negative and complex numbers.

The problem is has to do with the I π + Log[-1 + a] factor in the exponent. The other factors turn the solution into an expression containing a complex power of 1 = E^(2 I π).

FullSimplify[
 E^(((-1 + a^t) (I π + Log[-1 + a]))/(-1 + a)) // ComplexExpand,
 0 < a < 1 && k0 > 0 && t ∈ Integers && t > 0]

(* (1 - a)^((-1 + a^t)/(-1 + a)) E^((2 I (-1 + a^t) π)/(-1 + a)) *)

But it should be a real power of the positive base 1 - a:

FullSimplify[
 E^(I \[Pi] + Log[-1 + a]) // ComplexExpand, 
 0 < a < 1 && k0 > 0 && t \[Element] Integers && t > 0]

(* 1 - a *)

Almost anything you do to get rid of the extraneous factor seems like cheating, but the answer is erroneous anyway.

Simpler fix

My orginal answer was the result of spending time trying and failing to get Mathematica to look at k[t]^a as a real exponential instead of a complex exponential. But if we're going to cheat, then a simpler way is to replace the troublesome I π + Log[-1 + a] by Log[1 - a], which is what it should be:

rsol /. (I \[Pi] + Log[-1 + a]) -> Log[1 - a]

(* {k[t] -> (1 - a)^((-1 + a^t)/(-1 + a)) k0^a^t} *)

Original fix

Here's a way that gets and simplifies E^(I \[Pi] + Log[-1 + a]) within the expression using an extra TransformationFunctions:

Simplify[
 E^(((-1 + a^t) (I π + Log[-1 + a]))/(-1 + a)) k0^a^t,
 0 < a < 1 && k0 > 0 && t \[Element] Integers && t > 0,
 TransformationFunctions -> {
   Automatic,
   # /. Power[E, pow_] :> (
      ComplexExpand[E^#]^(pow/#) &@
        First @ Cases[pow, _?(! FreeQ[#, I] &)]) &}]

(* (1 - a)^((-1 + a^t)/(-1 + a)) k0^a^t *)
$\endgroup$
  • $\begingroup$ Thank you for your quick reply. Well, in the end I'll got to work with my manual solution. More worrying is that Mathematica's answer provides wrong numbers when we provide numeric values for the parameters. I am new to Mathematica, but that indeed seems like a bug. $\endgroup$ – benybenjamin Dec 23 '13 at 11:53
4
$\begingroup$

An alternative approach is to take the Log of both sides of the recurrence relation, solve the transformed relation, and exponentiate the result.

l[t] /. First@RSolve[{l[t + 1] == Log[1 - a] + l[t] a, l[0] == l0}, l[t], t];
FullSimplify[Exp[%] /. l0 -> Log[k0]]

(* (1 - a)^((-1 + a^t)/(-1 + a)) k0^a^t *)

By the way, the problem identified by MichaelE2 persists in version 11.0.1.

Workaround

A workaround based on a suggestion provided by Wolfram Technical Support for this particular problem is

SetOptions[PowerExpand, Assumptions -> 0 < a < 1];
FullSimplify[k[t] /. First@RSolve[{k[t + 1] == (1 - a) k[t]^a, k[0] == k0}, k[t], t],
    0 < a < 1]
SetOptions[PowerExpand, Assumptions -> Automatic]

(* (1 - a)^((-1 + a^t)/(-1 + a)) k0^a^t *)

PowerExpand, called by RSolve, needs to know that 0 < a < 1 to work properly here. Incidentally, SetOptions[PowerExpand, Assumptions -> True]; also seems to work.

$\endgroup$
1
$\begingroup$

The recursively defined function works if a numerical value of a (0

Here is a "cheat" using a wrapper then changing the form. It returns a form perhaps closer to the desired result:

k[n_, a_] := (1 - w[a]) k[n - 1, w@a]^w[a];
k[0, a_] := k0

then

FullSimplify[k[4, q] /. w -> ToExpression]

yields:

-(-1 + q) (-(-1 + q) (-(-1 + q) (-k0^q (-1 + q))^q)^q)^q

Continuing the cheat:

cheat[u_, var_] := Module[{a1, a2, a3, e1, e2},
  a1 = Simplify[u] /. (-1 + var) -> -h;
  a2 = a1 /. h -> 1;
  a3 = a1 /. k0 -> 1;
  e1 = Times @@ Cases[a2, Power[__, x_] :> x, {0, Infinity}];
  e2 = Times @@ Cases[a3, Power[_, x_] :> x, {0, Infinity}];
  (1 - var)^e2 k0^e1]

Tabulating from 1 to 6 (omitted 0 as case is defined and not accounted in cheat):

TableForm[
 Table[{j, cheat[k[j, m] /. w -> ToExpression, m]}, {j, 1, 6}]]

enter image description here

$\endgroup$
  • $\begingroup$ That works, thanks. However, it is a rather a cumbersome way for getting a solution for an easy problem. $\endgroup$ – benybenjamin Dec 23 '13 at 12:00
  • $\begingroup$ I agree entirely, hence my use of 'cheat' not really apt but am time poor at present...using Assuming and Assumptions did not work for me hence work around... $\endgroup$ – ubpdqn Dec 24 '13 at 0:21
  • $\begingroup$ When I defined recursively with condition on a I get correct form. Posted only to show symbolic manipulation could be forced. Look forward to insights from other answers. Happy holidays and peace to all. $\endgroup$ – ubpdqn Dec 24 '13 at 2:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.