Fix SymPy OverflowError in solow.md by using Rational

lectures/solow.md

@@ -3,10 +3,12 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.17.1
 kernelspec:
-  display_name: Python 3
-  language: python
   name: python3
+  display_name: Python 3 (ipykernel)
+  language: python
 ---
 
 (solow)=
@@ -117,16 +119,16 @@ The function $g$ from {eq}`solow` is then plotted, along with the 45-degree line
 Let's define the constants.
 
 ```{code-cell} ipython3
-A, s, alpha, delta = 2, 0.3, 0.3, 0.4
+A, s, α, δ = 2, 0.3, 0.3, 0.4
 x0 = 0.25
 xmin, xmax = 0, 3
 ```
 
 Now, we define the function $g$.
 
 ```{code-cell} ipython3
-def g(A, s, alpha, delta, k):
-    return A * s * k**alpha + (1 - delta) * k
+def g(A, s, α, δ, k):
+    return A * s * k**α + (1 - δ) * k
 ```
 
 Let's plot the 45-degree diagram of $g$.
@@ -139,7 +141,7 @@ def plot45(kstar=None):
 
     ax.set_xlim(xmin, xmax)
 
-    g_values = g(A, s, alpha, delta, xgrid)
+    g_values = g(A, s, α, δ, xgrid)
 
     ymin, ymax = np.min(g_values), np.max(g_values)
     ax.set_ylim(ymin, ymax)
@@ -202,11 +204,10 @@ If initial capital is above this level, then the reverse is true.
 Let's plot the 45-degree diagram to show the $k^*$ in the plot.
 
 ```{code-cell} ipython3
-kstar = ((s * A) / delta)**(1/(1 - alpha))
+kstar = ((s * A) / δ)**(1/(1 - α))
 plot45(kstar)
 ```
 
-
 From our graphical analysis, it appears that $(k_t)$ converges to $k^*$, regardless of initial capital
 $k_0$.
 
@@ -221,7 +222,7 @@ At this parameterization, $k^* \approx 1.78$.
 Let's define the constants and three distinct initial conditions
 
 ```{code-cell} ipython3
-A, s, alpha, delta = 2, 0.3, 0.3, 0.4
+A, s, α, δ = 2, 0.3, 0.3, 0.4
 x0 = np.array([.25, 1.25, 3.25])
 
 ts_length = 20
@@ -232,7 +233,7 @@ ymin, ymax = 0, 3.5
 ```{code-cell} ipython3
 def simulate_ts(x0_values, ts_length):
 
-    k_star = (s * A / delta)**(1/(1-alpha))
+    k_star = (s * A / δ)**(1/(1-α))
     fig, ax = plt.subplots(figsize=[11, 5])
     ax.set_xlim(xmin, xmax)
     ax.set_ylim(ymin, ymax)
@@ -243,7 +244,7 @@ def simulate_ts(x0_values, ts_length):
     for x_init in x0_values:
         ts[0] = x_init
         for t in range(1, ts_length):
-            ts[t] = g(A, s, alpha, delta, ts[t-1])
+            ts[t] = g(A, s, α, δ, ts[t-1])
         ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6,
                 label=r'$k_0=%g$' %x_init)
     ax.plot(np.arange(ts_length), np.full(ts_length,k_star),
@@ -319,14 +320,14 @@ levels of capital combine to yield global stability.
 To see this in a figure, let's define the constants
 
 ```{code-cell} ipython3
-A, s, alpha, delta = 2, 0.3, 0.3, 0.4
+A, s, α, δ = 2, 0.3, 0.3, 0.4
 ```
 
 Next we define the function $g$ for growth in continuous time
 
 ```{code-cell} ipython3
-def g_con(A, s, alpha, delta, k):
-    return A * s * k**alpha - delta * k
+def g_con(A, s, α, δ, k):
+    return A * s * k**α - δ * k
 ```
 
 ```{code-cell} ipython3
@@ -335,7 +336,7 @@ def plot_gcon(kstar=None):
     k_grid = np.linspace(0, 2.8, 10000)
 
     fig, ax = plt.subplots(figsize=[11, 5])
-    ax.plot(k_grid, g_con(A, s, alpha, delta, k_grid), label='$g(k)$')
+    ax.plot(k_grid, g_con(A, s, α, δ, k_grid), label='$g(k)$')
     ax.plot(k_grid, 0 * k_grid, label="$k'=0$")
 
     if kstar:
@@ -364,7 +365,7 @@ def plot_gcon(kstar=None):
 ```
 
 ```{code-cell} ipython3
-kstar = ((s * A) / delta)**(1/(1 - alpha))
+kstar = ((s * A) / δ)**(1/(1 - α))
 plot_gcon(kstar)
 ```
 
@@ -450,14 +451,14 @@ $$
 
 ```{code-cell} ipython3
 A = 2.0
-alpha = 0.3
-delta = 0.5
+α = 0.3
+δ = 0.5
 ```
 
 ```{code-cell} ipython3
 s_grid = np.linspace(0, 1, 1000)
-k_star = ((s_grid * A) / delta)**(1/(1 - alpha))
-c_star = (1 - s_grid) * A * k_star ** alpha
+k_star = ((s_grid * A) / δ)**(1/(1 - α))
+c_star = (1 - s_grid) * A * k_star ** α
 ```
 
 Let's find the value of $s$ that maximizes $c^*$ using [scipy.optimize.minimize_scalar](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize_scalar.html#scipy.optimize.minimize_scalar).
@@ -469,8 +470,8 @@ from scipy.optimize import minimize_scalar
 
 ```{code-cell} ipython3
 def calc_c_star(s):
-    k = ((s * A) / delta)**(1/(1 - alpha))
-    return - (1 - s) * A * k ** alpha
+    k = ((s * A) / δ)**(1/(1 - α))
+    return - (1 - s) * A * k ** α
 ```
 
 ```{code-cell} ipython3
@@ -510,21 +511,26 @@ plt.show()
 One can also try to solve this mathematically by differentiating $c^*(s)$ and solve for $\frac{d}{ds}c^*(s)=0$ using [sympy](https://www.sympy.org/en/index.html).
 
 ```{code-cell} ipython3
-from sympy import solve, Symbol
+from sympy import solve, Symbol, Rational
 ```
 
 ```{code-cell} ipython3
+# Use Rational for exact symbolic computation
+A = Rational(2)
+α = Rational(3, 10)
+δ = Rational(1, 2)
+
 s_symbol = Symbol('s', real=True)
-k = ((s_symbol * A) / delta)**(1/(1 - alpha))
-c = (1 - s_symbol) * A * k ** alpha
+k = ((s_symbol * A) / δ)**(1/(1 - α))
+c = (1 - s_symbol) * A * k ** α
 ```
 
 Let's differentiate $c$ and solve using [sympy.solve](https://docs.sympy.org/latest/modules/solvers/solvers.html#sympy.solvers.solvers.solve)
 
 ```{code-cell} ipython3
 # Solve using sympy
 s_star = solve(c.diff())[0]
-print(f"s_star = {s_star}")
+print(f"s_star = {float(s_star)}")
 ```
 
 Incidentally, the rate of savings which maximizes steady state level of per capita consumption is called the [Golden Rule savings rate](https://en.wikipedia.org/wiki/Golden_Rule_savings_rate).
@@ -576,23 +582,23 @@ Let's define the constants for lognormal distribution and initial values used fo
 
 ```{code-cell} ipython3
 # Define the constants
-sig = 0.2
-mu = np.log(2) - sig**2 / 2
+σ = 0.2
+μ = np.log(2) - σ**2 / 2
 A = 2.0
 s = 0.6
-alpha = 0.3
-delta = 0.5
+α = 0.3
+δ = 0.5
 x0 = [.25, 3.25] # list of initial values used for simulation
 ```
 
 Let's define the function *k_next* to find the next value of $k$
 
 ```{code-cell} ipython3
 def lgnorm():
-    return np.exp(mu + sig * np.random.randn())
+    return np.exp(μ + σ * np.random.randn())
 
-def k_next(s, alpha, delta, k):
-    return lgnorm() * s * k**alpha + (1 - delta) * k
+def k_next(s, α, δ, k):
+    return lgnorm() * s * k**α + (1 - δ) * k
 ```
 
 ```{code-cell} ipython3
@@ -604,7 +610,7 @@ def ts_plot(x_values, ts_length):
     for x_init in x_values:
         ts[0] = x_init
         for t in range(1, ts_length):
-            ts[t] = k_next(s, alpha, delta, ts[t-1])
+            ts[t] = k_next(s, α, δ, ts[t-1])
         ax.plot(np.arange(ts_length), ts, '-o', ms=4,
                 alpha=0.6, label=r'$k_0=%g$' %x_init)
 
@@ -621,7 +627,5 @@ def ts_plot(x_values, ts_length):
 ts_plot(x0, 50)
 ```
 
-
-
 ```{solution-end}
 ```