Difference Equations

Mathematica has a package that can help you solve difference equations. It is called RSolve (for recurrence solve). It has to be loaded into the kernel.

<< DiscreteMath`RSolve`

req  =  y[n + 1]  ==   -y[n]/2 + 6

RSolve[req, y[n], n]
TeXForm[
      (2^(-n)) (-4 (((-1))^n) + 
            (2^(2 + n)) + (((-1))^n) c)]

The equilibrium solution can be found by solving:

Solve[y == -y/2 + 6, y]

If we actually have an initial condition, we can solve the initial value difference equation in a manner similar to differential equations.

RSolve[{req, y[0] == 2}, y[n], n]

ysol[
        n_] = -2 (-2 + (((-(1/2)))^n))
data  =  Table[{n, ysol[n]} , {n, 0, 10}]
ListPlot[data, Frame -> True, 
      FrameLabel -> {"n", "y_n"}]
ListPlot[data, Frame -> True, 
      FrameLabel -> {"n", "y_n"}, 
      PlotJoined -> True]

If Mathematica can't solve it using RSolve, you can always generate the solution by iterating the difference equation. Note teh delayed assignments for y[0] and y[n].

y[0]  :=  2
y[n_]  :=   -y[n - 1]/2 + 6
data2  =  Table[{n, y[n]}, {n, 0, 10}]
ListPlot[data2]
TeXForm[data2]
TeXForm[N[data2]]

Plot[{x, -x/2 + 6}, {x, 0, 6}]

Discrete Logistic Equation, Computational Thoughts, and Chaos

Clear[y]

Before we go any further, we should note that we are not being very smart with our programming. By recursively defining things in the following manner, we are forcing Mathematica to climb down the tree to calculute each y[n]. For large n this is a lot of tree climbing! The ShowTime package will show us this.

<< Utilities`ShowTime`

Clear[y, data2]

y[0]  :=  1/3. 
y[n_]  :=   4 y[n - 1] ( 
          1 -  y[n - 1])
data2  =  Table[{n, y[n]}, {n, 0, 20}]
plot1  =  ListPlot[data2, PlotJoined -> True]

y[0]  :=  1/3. 
y[n_]  :=  y[n] =  
        4 y[n - 1] ( 1 -  y[n - 1])
data2  =  Table[{n, y[n]}, {n, 0, 20}]
plot1  =  ListPlot[data2, PlotJoined -> True]

The Discrete Logistic Equation

Clear[data2, y]
y[0]  :=  1/3. 
y[n_]  :=   3.65 y[n - 1] ( 
          1 -  y[n - 1])
data2  =  {{0, Evaluate[y[0]]}};
Do[{y[n]  =  Evaluate[y[n]], 
data2  = 
          Append[data2, {n, y[n]}]}, {n, 1, 40}]
plot1  =  ListPlot[data2, PlotJoined -> True]

Clear[data2, y]
y[0]  :=  1/3.   +  0.001
y[n_]  :=   3.65 y[n - 1] ( 
          1 -  y[n - 1])
data2  =  {{0, Evaluate[y[0]]}};
Do[{y[n]  =  Evaluate[y[n]], 
data2  = 
          Append[data2, {n, y[n]}]}, {n, 1, 40}]
plot2  =  
      ListPlot[data2, PlotJoined -> True, 
        PlotStyle -> {Hue[1], Dashing[{0.01, 0.01}]}]

Show[plot1, plot2, Frame -> True, 
      FrameLabel -> {"n", "y_n"}]

Clear[y]
RSolve[{y[n + 1] == 3.65 y[n] ( 1 -  y[n]), 
        y[0] == a}, y[n], n]

Here is a numerical solution to a differential equation.

Clear[f]
NDSolve[{f'[t] == f[t] ( 2.65 - 3.65 f[t]), 
        f[0] == 1/3}, f[t], {t, 0, 10}]
f[t_]  =  f[t]  /.  f[t]  ->  Part[%, 1, 1, 2]

Plot[f[t], {t, 0, 4}, PlotRange -> {0.3, 0.9}, Frame -> True, 
      FrameLabel -> {"t", "y(t)"}]

NonChaotic Regime

Clear[data2, y]
y[0]  :=  1/10
y[n_]  :=   2 y[n - 1] ( 1 -  y[n - 1])  +  
        y[n - 1]
data2  =  {{0, Evaluate[y[0]]}};
Do[{y[n]  =  Evaluate[y[n]], 
data2  = 
          Append[data2, {n, y[n]}]}, {n, 1, 10}]
N[data2]
plot1  =  
      ListPlot[data2, PlotJoined -> True, 
        PlotStyle -> {Hue[1], Dashing[{0.01, 0.01}]}]

Clear[f]
DSolve[{f'[t] == 2  f[t] ( 1 -  f[t]), 
        f[0] == 1/10}, f[t], t]

sol[
        t_]  = ((E^(2 t))/(9 + (E^(2 t))))
plot2  =  
      Plot[sol[t], {t, 0, 10}, PlotRange -> {0, 1.1}, 
        Frame -> True, 
        FrameLabel -> {"t", "y(t)"}]

Show[plot2, plot1, Frame -> True, 
      FrameLabel -> {"t", ""}]

Making the time delay smaller

Clear[data2, y]
y[0]  :=  1/10
y[n_]  :=   0.2 y[n - 1] ( 1 -  y[n - 1])  +  
        y[n - 1]
data2  =  {{0, Evaluate[y[0]]}};
Do[{y[n]  =  Evaluate[y[n]], 
data2  = 
          Append[data2, {n, y[n]}]}, {n, 1, 40}]
N[data2]
plot1  =  
      ListPlot[data2, PlotJoined -> True, 
        PlotStyle -> {Hue[1], Dashing[{0.01, 0.01}]}]

sol[
        t_]  = ((E^(0.2 t))/(9 + (E^(0.2 t))))
plot2  =  
      Plot[sol[t], {t, 0, 40}, PlotRange -> {0, 1.1}, 
        Frame -> True, 
        FrameLabel -> {"t", "y(t)"}]
Show[plot2, plot1, Frame -> True, 
      FrameLabel -> {"t", ""}]