plot(x * sin(1/x), x, -2, 2, plot_points=500)
def p(x, n):
return(taylor(sin(x), x, 0, n))
xmax = 15
n = 15
g = plot(sin(x), x, -xmax, xmax)
for d in range(n):
g += plot(p(x, 2 * (n-d) - 1), x, -xmax, xmax,
color=(1.7*d/(2*n), 1.5*d/(2*n), 1-3*d/(4*n)))
g.show(ymin=-2, ymax=2)
g.save('sine.png', ymin=-2, ymax=2)
f2(x) = 1
f1(x) = -1
f = piecewise( [ [(-pi,0),f1], [(0,pi),f2] ] )
S = f.fourier_series_partial_sum(20,pi)
g = plot(S, x, -8, 8, color='blue')
saw(x) = x - 2 * pi * floor((x + pi) / (2 * pi))
g += plot(saw(x) / abs(saw(x)), x, -8, 8, color='red')
g
t = var('t')
x = cos(t) + cos(7*t)/2 + sin(17*t)/3
y = sin(t) + sin(7*t)/2 + cos(17*t)/3
g = parametric_plot((x, y), (t, 0, 2*pi))
g.show(aspect_ratio=1)
t = var('t')
n = 20/19
g1 = polar_plot(1+2*cos(n*t),(t,0,n*36*pi),plot_points=5000)
g2 = polar_plot(1+1/3*cos(n*t),(t,0,n*36*pi),plot_points=5000)
g1.show(aspect_ratio=1)
g2.show(aspect_ratio=1)
Let us draw the curve given by the implicit equation $$\mathcal{C} = \left\{ z \in \mathbb{C} : |\cos(x^4)| = 1\right\}.$$
z = var('z')
g1 = complex_plot(abs(cos(z^4))-1, (-3,3), (-3,3), plot_points = 400)
g1.show(aspect_ratio=1)
f = lambda x, y: (abs(cos((x + I * y) ** 4)) - 1)
g2 = implicit_plot(f, (-3, 3), (-3, 3))
g2.show(aspect_ratio=1)
bar_chart([randrange(15) for i in range(20)])
bar_chart([x^2 for x in range(1,20)], width=0.2)
liste = [10 + floor(10*sin(i)) for i in range(100)]
bar_chart(liste)
finance.TimeSeries(liste).plot_histogram(bins=20)
Example. (Random walk) Starting from the origin O, a particle moves a distance $\ell$ every $t$ seconds, in a random direction, independently of the preceding moves. Let us draw an example of particle trajectory. The red line goes from the initial to the final position.
n, l, x, y = 10000, 1, 0, 0
p = [[0, 0]]
for k in range(n):
theta = (2 * pi * random()).n(digits=5)
x, y = x + l * cos(theta), y + l * sin(theta)
p.append([x, y])
g1 = line([p[n], [0, 0]], color='red', thickness=2)
g1 += line(p, thickness=.4)
g1.show(aspect_ratio=1)
Example. (Uniformly distributed sequences) Given a real sequence $(u_n)_{n\in \mathbb{N}^*}$, we construct the polygonal line whose successive vertices are the points in the complex plane $$z_N = \sum_{n\leq N} e^{2i \pi u_n}.$$
def UDS(u, length):
n = var('n')
z = lambda n: exp(2 * I * pi * u(n)).n()
vertices = [CC(0, 0)]
for n in range(1, length):
vertices.append(vertices[n - 1] + CC(z(n)))
line(vertices).show(aspect_ratio=1)
UDS(lambda n: n * sqrt(2), 200)
UDS(lambda n: n * ln(n) * sqrt(2), 10000)
UDS(lambda n: floor(n * ln(n)) * sqrt(2), 10000)
UDS(lambda n: nth_prime(n) * sqrt(2), 10000)
Example. (First-order linear differential equation) Let us draw the integral curves of the differential equation $$xy' - 2y = x^3.$$
# y' = (x^3 + 2y) / x
y = var('y')
vf = plot_vector_field((x, 2*y+x^3), (x,-2,2), (y,-1,1))
vf.show(ymin=-1, ymax=1)
x = var('x')
y = function('y')
DE = x*diff(y(x), x) == 2*y(x) + x^3
desolve(DE, [y(x),x])
# Symbolic solution.
x = var('x')
y = function('y')
sol = []
for i in [-2, -1.8..2]:
sol.append(desolve(DE, [y(x), x], ics=[1, i]))
sol.append(desolve(DE, [y(x), x], ics=[-1, i]))
g = plot(sol, x, -2, 2, color ='blue')
(g+vf).show(ymin=-1, ymax=1)
# Numerical solution.
x = var('x')
y = function('y')
DE = x*diff(y(x), x) == 2*y(x) + x^3
g = Graphics() # creates an empty graph
for i in [-2, -1.8..2]:
g += line(desolve_rk4(DE, y(x), ics=[ 1, i], step=0.05, end_points=[0, 2]))
g += line(desolve_rk4(DE, y(x), ics=[-1, i], step=0.05, end_points=[-2, 0]))
(g+vf).show(ymin=-1, ymax=1)
Example. (First-order non-linear differential equation) Let us draw the integral curves of the equation $$y'(t) + \cos(y(t)\cdot t) = 0.$$
y = var('y')
vf = plot_vector_field((1, -cos(x*y)), (x,-5,5), (y,-2,11))
vf.show()
import scipy
from scipy import integrate
f = lambda y, t: - cos(y * t)
t = [0, 0.1 .. 5]
p = Graphics()
for k in [0, 0.15 .. 10]:
y = integrate.odeint(f, k, t)
p += line(zip(t, y))
t = [0, -0.1 .. -5]
q = Graphics()
for k in [0, 0.15 .. 10]:
y = integrate.odeint(f, k, t)
q += line(zip(t, y))
g = p + q + vf
g.show()
Example. (Lotka-Volterra predator-prey model) We wish to represent graphically the variation of a set of prey and predators evolving according to a system of Lotka-Volterra equations: $$\begin{cases} \dfrac{du}{dt} &= au - buv\\ \dfrac{dv}{dt} &= -cv + dbuv\\ \end{cases}$$ where $u$ is the number of preys (for example rabbits), $v$ is the number of predators (for example foxes). In addition, the parameters $a$, $b$, $c$, $d$ describe the evolution of the populations: $a$ is the natural growth of rabbits without foxes to eat them, $b$ is the decrease of rabbits when foxes kill them, $c$ is the decrease of foxes without any rabbit to eat, and finally $d$ indicates how many rabbits are needed for a new fox to appear.
import scipy
from scipy import integrate
a, b, c, d = 1., 0.1, 1.5, 0.75
def dX_dt(X, t=0): # returns the population variation
return [a*X[0] - b*X[0]*X[1], -c*X[1] + d*b*X[0]*X[1]]
t = srange(0, 15, .01) # time scale
X0 = [10, 5] # initial conditions: 10 rabbits and 5 foxes
X = integrate.odeint(dX_dt, X0, t) # numerical solution
rabbits, foxes = X.T # shortcut for X.transpose()
p = line(zip(t, rabbits), color='red') # number of rabbits graph
p += text("Rabbits",(12,37), fontsize=10, color='red')
p += line(zip(t, foxes), color='blue') # idem for foxes
p += text("Foxes",(12,7), fontsize=10, color='blue')
p.axes_labels(["time", "population"])
p.show(gridlines=True)
n = 11
L = srange(6, 18, 12/n)
R = srange(3, 9, 6/n)
CI = list(zip(L, R)) # list of initial conditions
def g(x,y):
v = vector(dX_dt([x, y])) # for a nicer graph,
return v/v.norm() # we normalise the vector field
x, y = var('x, y')
q = plot_vector_field(g(x, y), (x, 0, 60), (y, 0, 36))
for j in range(n):
X = integrate.odeint(dX_dt, CI[j], t) # resolution
q += line(X, color=hue(.8-float(j)/(1.8*n))) # graph plot
q.axes_labels(["rabbits","foxes"])
q.show()
Example. (Evolute of the parabola) Let us find the equation of the evolute of the parabola $P$ of equation $y = x^2/4$, and show on the same graph the parabola $P$, some normals to $P$ and its evolute.
x, y, t = var('x, y, t')
alpha(t) = 1
beta(t) = t / 2
gamma(t) = t + t^3 / 8
env = solve([alpha(t) * x + beta(t) * y == gamma(t),
diff(alpha(t), t) * x + diff(beta(t), t) * y == diff(gamma(t), t),
], [x,y])
show(env)
f(x) = x^2 / 4
p = plot(f, -8, 8, rgbcolor=(0.2,0.2,0.4)) # the parabola
for u in [0, 0.1 .. 8]: # normals to the parabola
p += line([[u, f(u)], [-8*u, f(u) + 18]], thickness=.3)
p += line([[-u, f(u)], [8*u, f(u) + 18]], thickness=.3)
p += parametric_plot((env[0][0].rhs(),env[0][1].rhs()), \
(t, -8, 8), color='red') # draws the evolute
p.show(xmin=-8, xmax=8, ymin=-1, ymax=12, aspect_ratio=1)
t = var('t')
p = 2
x(t) = t
y(t) = t^2 / (2 * p)
f(t) = [x(t), y(t)]
df(t) = [x(t).diff(t), y(t).diff(t)]
d2f(t) = [x(t).diff(t, 2), y(t).diff(t, 2)]
T(t) = [df(t)[0] / df(t).norm(), df[1](t) / df(t).norm()]
N(t) = [-df(t)[1] / df(t).norm(), df[0](t) / df(t).norm()]
R(t) = (df(t).norm())^3 / (df(t)[0]*d2f(t)[1]-df(t)[1]*d2f(t)[0])
Omega(t) = [f(t)[0] + R(t)*N(t)[0], f(t)[1] + R(t)*N(t)[1]]
g = parametric_plot(f(t), (t,-8,8), color = 'green',thickness=2)
for u in srange(.4, 4, .2):
g += point(f(t=u), color='red', size = 40)
g += point(Omega(t=u), color='red', size = 10)
g += line([f(t=u), Omega(t=u)], color = 'brown', alpha = .5)
g += circle(Omega(t=u), R(t=u), color = 'blue')
g.show(aspect_ratio=1, xmin=-12, xmax=7, ymin=-3, ymax=12)
u, v = var('u, v')
h = lambda u, v: u^2 + 2*v^2
plot3d(h, (u,-1,1), (v,-1,1), aspect_ratio=[1,1,1])
Example. (A discontinuous function whose directional derivatives exist everywhere!) Study the existence in $(0, 0)$ of the directional derivatives and the continuity of the function $f$ from $\mathbb{R}^2$ to $\mathbb{R}$ defined by: $$f(x, y) = \begin{cases} \frac{x^2y}{x^4+y^2} & \text{if $(x,y) \neq (0,0)$ },\\ 0 & \text{if $(x,y) = (0,0)$ }. \end{cases}$$
f(x, y) = x^2 * y / (x^4 + y^2)
t, theta = var('t, theta')
limit(f( t*cos(theta), t*sin(theta) ) / t, t=0)
solve(f(x, y) == 1/2, y)
a = var('a')
h = f(x, a*x^2).simplify_rational()
h
plot(h, a, -4, 4)
p = plot3d(f(x,y),(x,-2,2),(y,-2,2),plot_points=[150,150])
for i in [-1 .. 1]:
p += plot3d( i / 4, (x, -2, 2), (y, -2, 2), \
color=hue( (i+2) / 10), opacity = .2)
p
Example. (Cassini surface) $$(a^2 + x^2 + y^2)^2 = 4a^2x^2 + z^4$$
x, y, z = var('x, y, z')
a = 1
h = lambda x, y, z:(a^2 + x^2 + y^2)^2 - 4*a^2*x^2-z^4
implicit_plot3d(h, (x,-3,3), (y,-3,3), (z,-2,2), plot_points=100)
Example. (A knot in space)
line3d([ ( -10*cos(t)-2*cos(5*t)+15*sin(2*t), \
-15*cos(2*t)+10*sin(t)-2*sin(5*t), \
10*cos(3*t), \
) for t in [0, 0.1 .. 6.4] ], radius=.5)