Jdi na navigaci předmětu

Konvergence číselné řady

var('k,n,x')

(k, n, x)

Sčítání řad

SageMath umí sečíst celou řadu konvergentních řad.

sum(2^(-k), k, 0, oo)

sum(1/(k*(k+1)), k, 1, oo)

show(sum(sin(k)/k, k, 1, oo).trig_reduce())

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, \pi - \frac{1}{2}

sum((-1)^k * 2^(2*k) / (factorial(2*k)), k, 0, oo)

cos(2)

V případě divergentní řady občas dokáže i "zakřičet".

sum(2^k, k, 0, oo)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
~/src/sage/local/lib/python3.8/site-packages/sage/interfaces/maxima_lib.py in sr_sum(self, *args)
    876         try:
--> 877             return max_to_sr(maxima_eval([[max_ratsimp],[[max_simplify_sum],([max_sum],[sr_to_max(SR(a)) for a in args])]]));
    878         except RuntimeError as error:

~/src/sage/local/lib/python3.8/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.EclObject.__call__ (build/cythonized/sage/libs/ecl.c:8599)()
    851         lispargs = EclObject(list(args))
--> 852         return ecl_wrap(ecl_safe_apply(self.obj,(<EclObject>lispargs).obj))
    853 

~/src/sage/local/lib/python3.8/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.ecl_safe_apply (build/cythonized/sage/libs/ecl.c:5897)()
    364     if error != NULL:
--> 365         raise RuntimeError("ECL says: {}".format(
    366             ecl_string_to_python(error)))

RuntimeError: ECL says: Error executing code in Maxima: sum: sum is divergent.

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-6-353789dc3396> in <module>
----> 1 sum(Integer(2)**k, k, Integer(0), oo)

~/src/sage/local/lib/python3.8/site-packages/sage/misc/functional.py in symbolic_sum(expression, *args, **kwds)
    572     """
    573     if hasattr(expression, 'sum'):
--> 574         return expression.sum(*args, **kwds)
    575     elif len(args) <= 1:
    576         return sum(expression, *args)

~/src/sage/local/lib/python3.8/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.sum (build/cythonized/sage/symbolic/expression.cpp:66042)()
  12407         """
  12408         from sage.calculus.calculus import symbolic_sum
> 12409         return symbolic_sum(self, *args, **kwds)
  12410 
  12411     def prod(self, *args, **kwds):

~/src/sage/local/lib/python3.8/site-packages/sage/calculus/calculus.py in symbolic_sum(expression, v, a, b, algorithm, hold)
    609 
    610     if algorithm == 'maxima':
--> 611         return maxima.sr_sum(expression,v,a,b)
    612 
    613     elif algorithm == 'mathematica':

~/src/sage/local/lib/python3.8/site-packages/sage/interfaces/maxima_lib.py in sr_sum(self, *args)
    882 # could not find an example where 'Pole encountered' occurred, though
    883 #            if "divergent" in s or 'Pole encountered' in s:
--> 884                 raise ValueError("Sum is divergent.")
    885             elif "Is" in s: # Maxima asked for a condition
    886                 self._missing_assumption(s)

ValueError: Sum is divergent.

Někdy ale může být výsledek okořeněný speciálními funkcemi.

sum(1/(k^2 + 1), k, 0, oo)

-1/2*I*psi(I) + 1/2*I*psi(-I)

psi?

Nebo si SageMath nemusí vědět rady, pak zadání nechá v symbolickém tvaru.

sum(1/(k^2 + k + 2^k), k, 0, oo)

sum(1/(k^2 + 2^k + k), k, 0, +Infinity)

Příklad konvergentní číselné řady

Na následujícím příkladě demonstrujeme definici konvergentní číselné řady

$\sum_{k=0}^\infty \frac{\sin(k)}{2^k}$

Pokud chceme rozhodnout o její konvergenci pomocí definice, musíme se pokusit najít explicitní výraz pro $n$ tý částečný součet,

$s_n = \sum_{k=0}^n \frac{\sin(k)}{2^k}.$

K tomu účelu využijeme SageMath.

s(n) = sum(sin(k)/2^k, k, 0, n)
show(s(n))

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{2^{n + 1} \sin\left(1\right) - \cos\left(n \arctan\left(\frac{\sin\left(1\right)}{\cos\left(1\right)}\right) + \arctan\left(\frac{\sin\left(1\right)}{\cos\left(1\right)}\right)\right) \sin\left(1\right) + {\left(\cos\left(1\right) - 2\right)} \sin\left(n \arctan\left(\frac{\sin\left(1\right)}{\cos\left(1\right)}\right) + \arctan\left(\frac{\sin\left(1\right)}{\cos\left(1\right)}\right)\right)}{2^{n} {\left(4 \, \cos\left(1\right) - 5\right)}}

V tomto výpočtu samozřejmě silně důvěřujeme Sage (resp. Maxima). Musíme její odpovědi proto prověřovat (následující ani předchozí kód samozřejmě není důkaz této součtové formule):

[
    (s(n) - sum([ sin(k)/2^k for k in range(n+1)])).simplify_full().reduce_trig()
    for n in range(5)
]

[0, 0, 0, 0, 0]

Nyní vypočteme limitu této posloupnosti částečných součtů $(s_n)_{n=0}^\infty$ , to už je jednoduché cvičení se znalostmi s BI-ZMA,

show(limit(s(n), n=oo))

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{2 \, \sin\left(1\right)}{4 \, \cos\left(1\right) - 5}

Jinak řečeno, naše řada je konvergentní a její součet je výše uvedené číslo,

$\sum_{k=0}^\infty \frac{\sin(k)}{2^k} = \frac{2\sin(1)}{5 - 4\cos(1)}$

Pojďme tuto situaci ilustrovat graficky.

list_plot(
    [ s(n) for n in range(21) ],
    gridlines=['automatic', [2*sin(1)/(5-4*cos(1))]],
    axes_labels=['$n$', '$s_n$']
)

Případně, chcete-li, podobně nepřesně numericky na deset cifer,

table([
    ['$n$', '$s_n\\approx$']] +
    [ [m, N(s(m), digits=10)] for m in range(21)
])

$n$	$s_n\approx$
$0$	$-2.563048209 \times 10^{-12}$
$1$	$0.4207354924$
$2$	$0.6480598491$
$3$	$0.6656998501$
$4$	$0.6183996942$
$5$	$0.5884333106$
$6$	$0.5840674434$
$7$	$0.5892001512$
$8$	$0.5930648319$
$9$	$0.5938697508$
$10$	$0.5933384802$
$11$	$0.5928502037$
$12$	$0.5927192045$
$13$	$0.5927704944$
$14$	$0.5928309563$
$15$	$0.5928508015$
$16$	$0.5928464084$
$17$	$0.5928390735$
$18$	$0.5928362087$
$19$	$0.5928364946$
$20$	$0.5928373653$

A přibližná hodnota skutečného součtu je, opět na deset cifer,

N(2*sin(1) / (5 - 4*cos(1)), digits=10)

0.5928376207

O této řadě lze samozřejmě velmi snadno ukázat, že je konvergentní srovnáním s geometrickou řadou. Nalézt součet je komplikovanější. Srovnání stojí na odhadu pro členy této řady

$\left| \frac{\sin(k)}{2^k} \right| \leq \frac{1}{2^k}, \quad k\in\mathbb{N}_0.$

Příklad divergentní číselné řady

Notoricky známou divergentní řadou je řada

$\sum_{k=1}^\infty \frac{1}{k}$

Vezměme její posloupnost částečných součtů

s(n) = sum(1/k, k, 1, n)

Její chování je znázorněno na následujícím grafu.

list_plot(
    [ s(n) for n in range(201) ],
    gridlines=True,
    axes_labels=['$n$', '$s_n$']
)

Zde pozorujeme kvalitativně jiné chování než v předchozím příkladě. Z přednášky máme i dokázáno, že $s_n \to +\infty$ .

Později během semestru ukážeme, že tato posloupnost do nekonečna ubíhá velmi pomalu, přesně jako přirozený logaritmus.

list_plot(
    [ s(n)/log(n) for n in range(2,201) ],
    gridlines=True, ymin=1/2, ymax=3/2,
    axes_labels=['$n$', '$s_n / \\ln(n)$']
)

Výpočty s využitím řad

Exponenciální funkce

Exponenciální funkci máme definovánu jako součet následující číselné řady

$e^x = \sum_{k=0}^\infty \frac{x^k}{k!}.$

O takto definované funkci jsme ukázali, že má všechny požadované vlastnosti. Navíc nám tato definice umožňuje přibližně počítat funkční hodnoty této funkce, například Eulerovo číslo, tedy $e^1$ .

s(n, x) = sum(x^k / factorial(k), k, 0, n)

Mimochodem,

s(0, 0).simplify()

0^0

Grafická ilustrace:

list_plot(
    [ s(n, 1).simplify() for n in range(21) ],
    gridlines=True,
    axes_labels=['$n$', '$s_n(1)$']
)

V tabulkové formě:

table([
    ['$n$', '$s_n(1)$', '$s_n(1)\\approx$']] + 
    [ [m, s(m, 1).simplify(), N(s(m, 1), digits=20)]
    for m in range(21)
])

$n$	$s_n(1)$	$s_n(1)\approx$
$0$	$1$	$1.0000000000000000000$
$1$	$2$	$2.0000000000000000000$
$2$	$\frac{5}{2}$	$2.5000000000000000000$
$3$	$\frac{8}{3}$	$2.6666666666666666667$
$4$	$\frac{65}{24}$	$2.7083333333333333333$
$5$	$\frac{163}{60}$	$2.7166666666666666667$
$6$	$\frac{1957}{720}$	$2.7180555555555555556$
$7$	$\frac{685}{252}$	$2.7182539682539682540$
$8$	$\frac{109601}{40320}$	$2.7182787698412698413$
$9$	$\frac{98641}{36288}$	$2.7182815255731922399$
$10$	$\frac{9864101}{3628800}$	$2.7182818011463844797$
$11$	$\frac{13563139}{4989600}$	$2.7182818261984928652$
$12$	$\frac{260412269}{95800320}$	$2.7182818282861685639$
$13$	$\frac{8463398743}{3113510400}$	$2.7182818284467590023$
$14$	$\frac{47395032961}{17435658240}$	$2.7182818284582297479$
$15$	$\frac{888656868019}{326918592000}$	$2.7182818284589944643$
$16$	$\frac{56874039553217}{20922789888000}$	$2.7182818284590422591$
$17$	$\frac{7437374403113}{2736057139200}$	$2.7182818284590450705$
$18$	$\frac{17403456103284421}{6402373705728000}$	$2.7182818284590452267$
$19$	$\frac{82666416490601}{30411275102208}$	$2.7182818284590452349$
$20$	$\frac{6613313319248080001}{2432902008176640000}$	$2.7182818284590452353$

Podle SageMath

N(e, digits=20)

2.7182818284590452354

Pozor, zamyslete se nad následujícím (také viz Mathematica notebook, zde je zajímavé oba programy porovnat):

N(exp(1.0), digits=20)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-27-0e6492cf627c> in <module>
----> 1 N(exp(RealNumber('1.0')), digits=Integer(20))

~/src/sage/local/lib/python3.8/site-packages/sage/misc/functional.py in numerical_approx(x, prec, digits, algorithm)
   1416         return numerical_approx_generic(x, prec)
   1417     else:
-> 1418         return n(prec, algorithm=algorithm)
   1419 
   1420 n = numerical_approx

~/src/sage/local/lib/python3.8/site-packages/sage/structure/element.pyx in sage.structure.element.Element.numerical_approx (build/cythonized/sage/structure/element.c:8182)()
    872         if prec is None:
    873             prec = digits_to_bits(digits)
--> 874         return numerical_approx_generic(self, prec)
    875 
    876     def n(self, prec=None, digits=None, algorithm=None):

~/src/sage/local/lib/python3.8/site-packages/sage/arith/numerical_approx.pyx in sage.arith.numerical_approx.numerical_approx_generic (build/cythonized/sage/arith/numerical_approx.c:2713)()
     63 
     64     if prec > inprec:
---> 65         raise TypeError("cannot approximate to a precision of %s bits, use at most %s bits" % (prec, inprec))
     66 
     67     # The issue is not precision, try conversion instead

TypeError: cannot approximate to a precision of 70 bits, use at most 53 bits

Výpočet $\pi$ (Bailey-Borwein-Plouffe formula, 1995)

Platí následující fascinující vztah

$\pi = \sum_{k=0}^\infty \frac{1}{16^k} \left( \frac{4}{8k+1} - \frac{2}{8k+4} - \frac{1}{8k+5} - \frac{1}{8k+6} \right).$

s(n) = sum(1/16^k * (4 / (8*k+1) - 2 / (8*k+4) - 1 / (8*k+5) - 1 / (8*k+6)), k, 0, n)
show(s(n))

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, {\sum_{k=0}^{n} \frac{2 \, {\left(120 \, k^{2} + 151 \, k + 47\right)}}{{\left(512 \, k^{4} + 1024 \, k^{3} + 712 \, k^{2} + 194 \, k + 15\right)} 16^{k}}}

list_plot(
    [ s(n).simplify() for n in range(21) ],
    gridlines=True,
    axes_labels=['$n$', '$s_n$']
)

table([
    ['$n$', '$s_n$', '$s_n\\approx$']] +
    [ [m, s(m).simplify(), N(s(m), digits=30)]
    for m in range(21)
])

$n$	$s_n$	$s_n\approx$
$0$	$\frac{47}{15}$	$3.13333333333333333333333333333$
$1$	$\frac{102913}{32760}$	$3.14142246642246642246642246642$
$2$	$\frac{615863723}{196035840}$	$3.14158739034658152305211128741$
$3$	$\frac{357201535487}{113700787200}$	$3.14159245756743538183700455506$
$4$	$\frac{16071212445820879}{5115625817702400}$	$3.14159264546033631955702122244$
$5$	$\frac{40413742330349316707}{12864093722915635200}$	$3.14159265322808753473437803554$
$6$	$\frac{4318127540987083098959311}{1374502686106089789849600}$	$3.14159265357288082778524076190$
$7$	$\frac{16331158360096799798177512637}{5198369158853231585211187200}$	$3.14159265358897270494077776717$
$8$	$\frac{2090388270092909093421859664191}{665391252333213642907031961600}$	$3.14159265358975227523617786840$
$9$	$\frac{64251934196540737654784844866951}{20452025861189303550405613977600}$	$3.14159265358979114638877696591$
$10$	$\frac{265231984363320332474665054275863923}{84425962754989445056074374499532800}$	$3.14159265358979312961417056404$
$11$	$\frac{4437861562367075948602657783785247810411}{1412615208816483394678236434126182809600}$	$3.14159265358979323271129226193$
$12$	$\frac{2782574702496655561202854740689196451804989707}{885721036849605620330411670088589631081676800}$	$3.14159265358979323815476632250$
$13$	$\frac{2426405140577083649593805922597689130605532212899}{772348744132856100928118976317250158303222169600}$	$3.14159265358979323844597750194$
$14$	$\frac{517658978311277334141536972885988480280412713877204531}{164775970468280051996408614883427417773042630551142400}$	$3.14159265358979323846173248204$
$15$	$\frac{28471243807120253377792333689239912281715167328354127149}{9062678375755402859802473818588507977517344680312832000}$	$3.14159265358979323846259317467$
$16$	$\frac{9278436701814805132782128093787327904238201772972845154896687}{2953418130518176727175308587792171687777172423177788194816000}$	$3.14159265358979323846264059514$
$17$	$\frac{722010830388420876212574684706667202761999027022892733523938606167}{229823185244402440201873813067635632056068449282002766167801856000}$	$3.14159265358979323846264322742$
$18$	$\frac{3442547639291990737781556257863191176805339146269075536550349546566411}{1095796947245310834882534340706486693643334366176589189088079249408000}$	$3.14159265358979323846264337452$
$19$	$\frac{13137821036965711114062966885739518573991633324975949215892604857282823589}{4181898318979566241553253362394632141925623428827149449925983670894592000}$	$3.14159265358979323846264338278$
$20$	$\frac{69788105348361857437902480107428396870369348283327260905420704007804969465293}{22214243870419455875130881861040285937908911653929817878006825259792072704000}$	$3.14159265358979323846264338325$

N(pi, digits=30)

3.14159265358979323846264338328