Zaporedje je funkcija iz naravnih števil v neko množico $\mathcal{Z}$.
$$ a:\mathbb{N}\to\mathcal{Z}$$
In [1]:
# zaporedje definiramo kot funkcijo
a = lambda n: n**10/2**n
for n in range(10):
print("%f" % a(n))
In [2]:
from matplotlib import pyplot as plt
%matplotlib inline
n = range(30)
plt.plot(n,[a(k) for k in n],'*')
#plt.semilogy(n,[a(k) for k in n],'*')
plt.title("prvih %d členov zaporedja" % len(n))
plt.show()
Zaporedje je padajoče, če je naslednji člen manjši od prejšnjega: $a_n>a_{n+1}$. Iz grafa vidimo, da je od približno 15 člena naprej zaporedje padajoče. Poiščimo $n$ pri katerem se to zgodi še s programom.
In [3]:
n = 0
while a(n)<a(n+1):
n += 1
n
Out[3]:
Graf in program seveda nista dokaz, da je zaporedje padajoče za vse $n>14$. Lahko pa to preverimo, tako da rešimo neenačbo $$ a_n>a_{n+1}$$ $$ \frac{n^{10}}{2^n}>\frac{(n+1)^{10}}{2^{n+1}}$$
In [4]:
import sympy as sym
sym.init_printing() # lepši izpis formul
from sympy.solvers.inequalities import solve_univariate_inequality
from sympy import Symbol
n = Symbol('n', real=True) # n simbol, ki je realno število
neenacba = a(n) > a(n+1)
sym.simplify(neenacba)
Out[4]:
In [6]:
solve_univariate_inequality(neenacba,n)
Neenačba je previsoke stopnje, zato je SymPy ne zna rešiti. Rešitev lahko poiščemo sami, tako da najprej rešimo enačbo. Nato določimo, na katerih intervalih neeančba velja.
In [7]:
from sympy.solvers import solve
resitve = solve(neenacba.lhs-neenacba.rhs)
resitve
Out[7]:
In [8]:
[float(x) for x in resitve]
Out[8]:
In [9]:
neenacba.subs(n,15) # v enačbo vstavimo n=15
Out[9]:
In [10]:
# narišemo nekaj prvih členov
n0 = 20
a = n0*[0]
a[0] = 1/2
for n in range(n0-1):
a[n+1] = 1+1./a[n]
print("Prvih %d členov zaporedja\n" % n0)
print((n0*"%0.8f\n")%tuple(a))
In [11]:
from matplotlib import pyplot as plt
%matplotlib inline
plt.plot(range(n0),a,'o')
plt.title("Členi zaporedja $a_n$")
plt.show()
Prvih nekaj členov nam razkrije, da je zaporedje lihih členov $a_{2n+1}$ padajoče, zaporedje sodih členov $a_{2n}$ pa naraščajoče. Da bi to preverili/dokazali, zapišimo rekurzivno enačbo, ki ji zadoščata podzaporedji $a_{2n+1}$ in $a_{2n}$.
Rekurzivna enačba
$$a_{n+1}=f(a_n)=1+\frac{1}{a_n}$$je določena s funkcijo $f(x)=1+\frac{1}{x}$. Za zaporedji sodih in lihih členov $a_{2n}$ in $a_{2n+1}$ dobimo rekurzivni formuli, če dvakrat uporabimo funkcijo $f$
$$a_{2(n+1)}=f(a_{2n+1})=f(f(a_{2n}))$$in
$$a_{2(n+1)+1} = f(f(a_{2n+1}))$$
In [12]:
import sympy as sym
sym.init_printing()
# zaporedje je podano s funkcijo f(x) = 1+1/x
f = lambda x: 1+1/x
n = sym.Symbol('n')
a_n = sym.Symbol('a_n')
f(f(a_n))
Out[12]:
In [13]:
sym.simplify(f(f(a_n)))
Out[13]:
Zaporedji lihih členov in sodih členov zadoščata rekurzivni formuli $$a_{n+2}=\frac{2 a_{n} + 1}{a_{n} + 1}.$$ Zaporedje $a_{2n}$ bo naraščalo, če bo veljalo $$a_{2n+2} > a_{2n}$$ oziroma $$\frac{2 a_{2n} + 1}{a_{2n} + 1} > a_{2n}.$$
In [14]:
ff = lambda x: f(f(x))
a_2n = sym.Symbol('a_2n')
sym.simplify(ff(a_2n) > a_2n)
Out[14]:
In [15]:
sym.solve_univariate_inequality(ff(a_2n)>a_2n,a_2n)
Out[15]:
Zaporedje sodih členov $a_{2n}$ bo padajoče, če so členi $$a_{2n}>\frac{1}{2}+\frac{\sqrt{5}}{2}$$ in naraščajoče, če so členi $$a_{2n}<\frac{1}{2}+\frac{\sqrt{5}}{2}.$$ Podobno velja za zaporedje lihih členov $a_{2n+1}$.
In [16]:
neenacba = ff(a_2n)<1/sym.Integer(2)+sym.sqrt(sym.Integer(5))/2
neenacba
Out[16]:
In [17]:
sym.solve_univariate_inequality(neenacba,a_2n)
Out[17]:
Vidimo, da je $$a_{2n}<\frac{1}{2}+\frac{\sqrt{5}}{2}\implies a_{2n+2}<\frac{1}{2}+\frac{\sqrt{5}}{2}.$$ Po indukciji sklepamo, da so vsi členi $a_{2n}<\frac{1}{2}+\frac{\sqrt{5}}{2}$, če to velja za prvi člen $$a_0= \frac{1}{2}<\frac{1}{2}+\frac{\sqrt{5}}{2}.$$ Ker so vsi členi $a_{2n}<\frac{1}{2}+\frac{\sqrt{5}}{2}$, je zaporedje sodih členov naraščajoče. Podobno dokažemo, da je zaporedje lihih členov $a_{2n+1}$ padajoče.
In [18]:
resitve = sym.solve(sym.Eq(a_n,f(a_n)),a_n)
resitve
Out[18]:
Rešitvi sta nam že poznani iz predhodne analize. Limita je očitno enaka prvi rešitvi $$\frac{1}{2}+\frac{\sqrt{5}}{2},$$ saj je $a_n>0$.
In [19]:
resitve[0].evalf() # numerični približek za limito se ujema s prej izračunanimi členi zaporedja
Out[19]:
Rekurzivno podano zaporedje z rekurzivno formulo $$a_{n+1}=f(a_n)$$ si lahko predstavljamo tudi grafično. Limita zaporedja $a_n$ je rešitev enačbe $$a=f(a)$$ Grafično prestavlja limita presečišče grafa funkcije $y=f(x)$ in simetrale lihih kvadrantov $y=x$.
In [21]:
import numpy as np
t=np.linspace(0.4,4)
l1, = plt.plot(t,t,label="$y=x$")
l2, = plt.plot(t,f(t),label="$y=f(x)$")
x = sym.Symbol('x')
resitve = sym.solve(sym.Eq(x,f(x)),x) # limita/rešitev enačbe x=f(x)
plt.plot(resitve[0],f(resitve[0]),'o')
# oznake in legenda
plt.annotate("$x=f(x)$",xy=(resitve[0],resitve[0]),xytext=(resitve[0]+0.3,resitve[0]))
plt.title("Fiksna točka funkcije $f(x)$")
legend1 = plt.legend(handles=[l1],loc=3)
plt.gca().add_artist(legend1)
plt.legend(handles=[l2],loc=2)
plt.axes().set_aspect('equal', 'datalim')
plt.show()
In [22]:
# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation.JSAnimation import IPython_display
import matplotlib.animation as animation
t=np.linspace(0.4,4)
fig = plt.figure()
plt.plot(t,t)
plt.plot(t,f(t))
plt.title("Konstrukcija zaporedja $a_{n+1}=f(a_n)$")
plt.ylim(-0.05,4)
plt.axes().set_aspect('equal', 'datalim')
ims = []
pts = []
sled = []
for i in range(4):
if a[i]>resitve[0]:
dx1,dx2 = -0.2,1
else:
dx1,dx2 = 1,-0.2
pts = pts + plt.plot(a[i],0,'ob') + [plt.annotate('$a_%d$'%i,xy=(a[i], 0), xytext=(a[i], 0.2))]
ims.append(sled+pts)
vl = plt.plot([a[i],a[i]],[0,f(a[i])],'-ro')
vl = vl + [plt.annotate('$(a_%d,f(a_%d))$'%(i,i),xy=(a[i], f(a[i])), xytext=(a[i]-dx1, f(a[i])+0.2))]
ims.append(sled+pts+vl)
hl = plt.plot([a[i],f(a[i])],[f(a[i]),f(a[i])],'-ro')
hl = hl + [plt.annotate('$(a_%d=f(a_%d))$'%(i+1,i),xy=(f(a[i]), f(a[i])), xytext=(f(a[i])-dx2, f(a[i])+0.2))]
ims.append(sled+pts+vl+hl)
#ims.append(hl)
vll = plt.plot([f(a[i]),f(a[i])],[f(a[i]),0],'-ro')
sled = sled + plt.plot([a[i],a[i],f(a[i]),f(a[i])],[a[i],f(a[i]),f(a[i]),f(f(a[i]))],'-',color='0.75')
ims.append(sled+pts+vl+hl+vll)
pts = pts + plt.plot(a[i+1],0,'ob') + [plt.annotate('$a_%d$'%i,xy=(a[i+1], 0), xytext=(a[i+1], 0.2))]
ims.append(pts)
animation.ArtistAnimation(fig,ims,interval=800,blit=True,repeat_delay=3000)
Out[22]:
Zaporedje grafično konstruiramo v več korakih. Za vsak nov člen zaporedja $a_{n+1}$ izvedemo naslednje korake:
Zaprtje (ang. closure) omogoča, da funkciji dodamo „spomin“, ne da bi pri tem uporabili globalne spremenljivke. Funkcijo, ki računa rekurzivno zaporedje, bomo popravili, da si bo zapomnila že izračunane člene.
In [ ]:
def gen_rekurzijo(a0,fun):
"""Funkcija definira rekurzivno zaporedje z rekurzivno formulo
a(n+1)=fun(a(n)) in začetnim členom
a(0)=a0"""
zaporedje = [a0]
def a(n):
if len(zaporedje)-1<n: # če člen zaporedja še ni izračunam
vrednost = fun(a(n-1)) # rekurzivna formula
zaporedje.append(vrednost) # shranimo vrednost za kasneje
return vrednost
else:
print("Že izračunani členi:", zaporedje)
return zaporedje[n] # člen a(n) je že izračunan
return a
In [ ]:
koren2 = gen_rekurzijo(2,lambda x:(x+2/x)/2)
koren2(1)
koren2(3)
In [ ]:
n = 0
while abs(koren2(n)**2-2)>1e-12:
n += 1
print("n: %d, koren: %1.12f" % (n,koren2(n)))
In [ ]:
import disqus
%reload_ext disqus
%disqus matpy
In [ ]: