Graphing

May 26, 2021

[1]: import numpy as np

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inline

[9]: x = np.linspace(0, 2*np.pi, 10000)

f1 = 2*np.sin(x)
f2 = 2*np.sin(x - np.pi/18)
f3 = 2*np.sin(3*x - np.pi/18)

fig, ax = plt.subplots(figsize=(12,6))

plt.plot(x, f1, label="$2 sin(x)$")

plt.plot(x, f2, label="$2 sin(x - \pi / 18)$")
plt.plot(x, f3, label="$2 sin(3(x - \pi / 18))$")

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

plt.legend()
plt.grid()

plt.show()

1
[40]: x = np.linspace(-15, 15, 10000)

f1 = 1000 - x**2
f2 = -20*x+1100

fig, ax = plt.subplots(figsize=(12,6))

plt.plot(x, f1, label="$f_1$")

plt.plot(x, f2, label="$f_2$")
plt.vlines(0, 1000, 1100, colors='k', linestyle='dashed')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

plt.legend()
plt.grid()

plt.show()

2
0.0.1 Example - Power law model

The table shows the mean (average) distances d of the planets from the sun (taking the unit of
measurement to be the distance from the earth to the sun) and their periods T (time of revolution
in years)
Planet
d
T
Mercury
0.387
0.241
Venus
0.723
0.615
Earth
1.000
1.000
Mars
1.523
1.881

3
 Jupiter
5.203
11.861
Saturn
9.541
29.457
Uranus
19.190
84.008
Neptune
30.086
164.784
[2]: d = np.array([0.387, 0.723, 1, 1.523, 5.203, 9.541, 19.19, 30.086])
T = np.array([0.241, 0.615, 1, 1.881, 11.861, 29.457, 84.008, 164.784])

x = np.linspace(0, 35, 100)

def func(x, a, n):

return a*x**n

param, _ = curve_fit(func, d, T)

textstr = '$T(d)$ = %.3f$x^{%.3f}$'%(param[0], param[1])

fig, ax = plt.subplots(figsize=(12,6))
plt.scatter(d, T, marker='o', facecolors='none', edgecolors='b', s=80,␣
,→label="data")

plt.plot(x, func(x, *param), '--', color='r', label ="fit")

plt.ylabel("$T$")
plt.xlabel("$d$")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.legend()
plt.legend(bbox_to_anchor=(1,0.5), loc="center left", frameon=False)

ax.text(0.25, 0.65, textstr, transform=ax.transAxes, fontsize=14,

verticalalignment='top')

plt.show()

4
 0.0.2 Example

What happens to the graph of the equation y 2 = cx3 + x2 as c varies?

[3]: cs = [1, 0.5, 0.1, -1]
domains = [np.linspace(-5, -1/c, 100) if c < 0 else np.linspace(-1/c, 5, 100)␣
,→for c in cs]

cmap = plt.get_cmap('gnuplot')
colors = [cmap(i) for i in np.linspace(0, 1, len(cs)+1)]

fig, ax = plt.subplots(figsize=(12,6))

for i, c in enumerate(cs):
x = domains[i]
plt.plot(x, np.sqrt(c*x**3+x**2), color=colors[i], label ="c=%.2f"%cs[i])
plt.plot(x, -np.sqrt(c*x**3+x**2), color=colors[i])

textstr = '$y^2=cx^3+x^2$'
ax.text(0.05, 0.85, textstr, transform=ax.transAxes, fontsize=14,
verticalalignment='top')

plt.legend()
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.legend(bbox_to_anchor=(1,0.6), loc="upper left", frameon=False)

plt.show()

5
 0.0.3 Example

Estimate the values of x such that ex > 1, 000, 000, 000

[6]: x = np.linspace(0, 21, 1000)
y = np.exp(x)
y_base = np.array([1000000000 for _ in x])

fig, ax = plt.subplots(figsize=(12,6))

plt.plot(x, y)
plt.plot(x, y_base)

idx = np.argwhere(np.diff(np.sign(y - y_base))).flatten()

plt.vlines(x[idx], -1, 1.3e9, colors='k', linestyle='dashed')

textstr = '$x=%.3f$'%x[idx]
ax.text(0.97, 0.76, textstr, transform=ax.transAxes, fontsize=14,
verticalalignment='top')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

plt.show()

6
 0.0.4 Example

Under ideal conditions a certain bacteria population is known to double every three hours. Suppose
that there are initially 100 bacteria.
[5]: p0 = 100
t = np.linspace(0, 28, 100)

# population at time t
p = 2**(t/3) * p0

# population at t = 15 hours
p15 = [2**(15/3) * p0 for _ in t]

# population at t = 20 hours
p20 = [2**(20/3) * p0 for _ in t]

p_base = [50000 for _ in t]

idx = np.argwhere(np.diff(np.sign(p - p_base))).flatten()
# population at t = 26
p26 = [2**(t[idx]/3) * p0 for _ in t]

fig, ax = plt.subplots(figsize=(12,6))

plt.plot(t, p)
plt.plot(t, p15, 'k--')
plt.plot(t, p20, 'k--')
plt.plot(t, p26, 'k--')

7
textstr1 = 'Population after 20 hours = %d'%(100*2**(20/3))
ax.text(0.97, 0.2, textstr1, transform=ax.transAxes, fontsize=10,
verticalalignment='top')

textstr2 = 'Population after 15 hours = %d'%(100*2**(15/3))

ax.text(0.97, 0.1, textstr2, transform=ax.transAxes, fontsize=10,
verticalalignment='top')

textstr3 = 'Population after %.2f hours = %d'%(t[idx], (100*2**(t[idx]/3)))

ax.text(0.97, 0.76, textstr3, transform=ax.transAxes, fontsize=10,
verticalalignment='top')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

plt.xlabel("Time (hours)")
plt.ylabel("Population")
plt.show()

0.0.5 Example

The table gives the population of Malaysia, in millions, for the years 1950-2000. Estimate the
population in 1975 and predict the population in the years 2010 and 2020
Year
Population
1950
6.1

8
 1955
7.0
1960
8.1
1965
9.5
1970
10.9
1975
12.3
1980
13.8
1985
15.7
1990
17.8
1995
20.4
2000
23.0
[4]: y0 = 1950
y = np.linspace(1950, 2000, 11)
p = np.array([6.1, 7.0, 8.1, 9.5, 10.9, 12.3, 13.8, 15.7, 17.8, 20.4, 23.0])

x = np.linspace(1940, 2025, 100)

def func(x, a, e, c, d):

return a*e**(c*(y0 - x)) + d

param, _ = curve_fit(func, y, p)

# population in 1975
p1975 = [func(1975, *param) for _ in x]

# population 2010
p2010 = [func(2010, *param) for _ in x]

9
# population 2020
p2020 = [func(2020, *param) for _ in x]

fig, ax = plt.subplots(figsize=(12,6))

plt.scatter(y, p, marker='o', facecolors='none', edgecolors='b', s=80,␣

,→label="data")

plt.plot(x, func(x, *param), '--', color='r', label ="fit")

plt.plot(x, p1975, 'k--')
plt.plot(x, p2010, 'k--')
plt.plot(x, p2020, 'k--')

textstr1 = 'Population in 1975 = %.1f (actual 12.3)'%(func(1975, *param))

ax.text(0.97, 0.25, textstr1, transform=ax.transAxes, fontsize=10,
verticalalignment='top')

textstr2 = 'Population in 2010 = %.1f (actual: 28.2)'%(func(2010, *param))

ax.text(0.97, 0.66, textstr2, transform=ax.transAxes, fontsize=10,
verticalalignment='top')

textstr3 = 'Population in 2020 = %.1f (actual 32.4)'%(func(2020, *param))

ax.text(0.97, 0.85, textstr3, transform=ax.transAxes, fontsize=10,
verticalalignment='top')

plt.xlabel("Year")
plt.ylabel("Population (in millions)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

plt.legend()
plt.legend(bbox_to_anchor=(1,0.45), loc="center left", frameon=False)

textstr = '$p(y)$ = %.3f*$%.3f^{%.3f (%d - y)} %.3f$'%(param[0], param[1],␣

,→param[2], y0, param[3])

ax.text(0.05, 0.45, textstr, transform=ax.transAxes, fontsize=14,

verticalalignment='top')

plt.show()

10
[ ]:

11