---

Κεφάλαιο 3

---

3. Λογισμός συναρτήσεων

3.1 Ορισμός μαθηματικών συναρτήσεων στο Sage

Στην Python μια συνάρτηση ορίζεται με την εντολή def. Χρησιμοποιώντας τον ορισμό αυτό μπορούμε να υπολογίζουμε τις τιμές της συνάρτησης σε διάφορα σημεία, καθώς επίσης και να σχεδιάζουμε την γραφική της παράσταση. Με την βοήθεια του Sage μπορούμε επιπλέον να παραγωγίσουμε, να ολοκληρώσουμε, να βρούμε όρια της συνάρτησης, να αναπτύξουμε την συνάρτηση σε σειρά Taylor γύρω από ένα σημείο, και γενικότερα να εφαρμόσουμε ένα πλήθος από μαθηματικές έννοιες για την μελέτη της, καθαρά συμβολικά!
Επίσης, το Sage μας παρέχει έναν πιο άμεσο τρόπο ορισμού μιας συνάρτησης

f(x,y) = sin(x+y) + cos(x-y)
print f, type(f)

(x, y) |--> cos(x - y) + sin(x + y) <type 'sage.symbolic.expression.Expression'>

Παρατηρήστε τον τρόπο που το Sage εμφανίζει την συνάρτηση, καθώς και ότι η συνάρτηση $f(x,y)$ είναι κι αυτή μια συμβολική έκφραση. Ας υπολογίσουμε μια τιμή της συνάρτησης με το Sage, κι ας ζητήσουμε μετά το αόριστο ολοκλήρωμά της ως προς $x$.

print f(3*pi/4,pi)
sage: print f.integrate(x)

0
(x, y) |--> -cos(x + y) - sin(-x + y)

Παρατηρούμε ότι το αποτέλεσμα της ολοκλήρωσης της $f(x,y)$ ως προς την μεταβλητή $x$, εμφανίζεται και πάλι με την μορφή μιας συνάρτησης.
Ας δούμε κάποιες διαφορές μεταξύ του Sage και της Python, όσον αφορά στον ορισμό συναρτήσεων

def g(x,y):
    return sin(x+y) + cos(x-y)
print g, type(g)

<function g at 0x7f231f2656e0> <type 'function'>

Όμως δεν μπορούμε να ολοκληρώσουμε την g(x,y) γιατί όπως μας πληροφορεί το Sage, to αντικείμενο 'function' δεν έχει την ιδιότητα 'integrate'

g.integrate(x)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-4-55863eb54981> in <module>()
----> 1 g.integrate(x)

AttributeError: 'function' object has no attribute 'integrate'

Επειδή όμως στο Sage έχουμε την δυνατότητα να ορίζουμε συναρτήσεις, και το Sage βλέπει τα πάντα συμβολικά, μπορούμε να θέσουμε την τιμή της g σε μια συμβολική μεταβλητή e

e = g(x,y)
print e, type(e)

cos(x - y) + sin(x + y) <type 'sage.symbolic.expression.Expression'>

Τώρα μπορούμε να μετατρέψουμε την συμβολική μεταβλητή e σε μια συνάρτηση, ας την πούμε h(x,y)

h(x,y) = e
print h, type(h)
print h.integrate(x)

(x, y) |--> cos(x - y) + sin(x + y) <type 'sage.symbolic.expression.Expression'>
(x, y) |--> -cos(x + y) - sin(-x + y)

απ' όπου παρατηρούμε ότι η h(x,y) είναι ορισμένη ως συμβολική έκφραση, και συνεπώς μπορούμε να την ολοκληρώσουμε.

3.2 Αλματικές συναρτήσεις

Η απλούστερη τμηματικά ασυνεχής συνάρτηση είναι η αλματική, δηλαδή η συνάρτηση εκείνη που η τιμή της αλλάζει από 0 σε 1, στο σημείο $x=0$. Ισοδύναμα λέμε ότι παρουσιάζει ένα άλμα ύψους μονάδας στο σημείο $x=0$. Υπάρχει άραγε η συνάρτηση αυτή στο Sage;

print unit_step, type(unit_step)
print 'at -3 :', unit_step(-3)
print 'at 4 :', unit_step(4)

unit_step <class 'sage.functions.generalized.FunctionUnitStep'>
at -3 : 0
at 4 : 1

Η συνάρτηση unit_step είναι ορισμένη στο Sage στην κλάση sage.functions.generalized.FunctionUnitStep.
Μπορούμε να δούμε τη βοήθεια που μας προσφέρει το Sage για την συνάρτηση αυτή με την εντολή

print sage.functions.generalized.FunctionUnitStep.__doc__

    The unit step function, `\mathrm{u}(x)` (``unit_step(x)``).

    INPUT:

    -  ``x`` - a real number or a symbolic expression

    DEFINITION:

    The unit step function, `\mathrm{u}(x)` is defined in Sage as:

        `\mathrm{u}(x) = 0` for `x < 0` and `\mathrm{u}(x) = 1` for `x \geq 0`

    EXAMPLES::

        sage: unit_step(-1)
        0
        sage: unit_step(1)
        1
        sage: unit_step(0)
        1
        sage: unit_step(x)
        unit_step(x)
    

Συνεπώς η συνάρτηση unit_step στο Sage ορίζεται ως εξής $$ $$ $$ \mathrm{u}(x) = \left\lbrace \begin{array}{cc} 1\,, & x \geq 0 \,, \\ 0\, & x < 0\,. \end{array} \right.$$

p = plot(unit_step, -1, 1, exclude=[0],thickness=2)
p.show(figsize=3)

Το αξιοσημείωτο είναι ότι το Sage αναγνωρίζει την αλματική συνάρτηση στην κλάση των γενικευμένων συναρτήσεων (ή αλλιώς κατανομών), οι οποίες με μια γενικευμένη έννοια, έχουν παραγώγους οποιασδήποτε τάξης! Για παράδειγμα, η πρώτη παράγωγος της "συνάρτησης" unit_step, ή αλλιώς της κατανομής Heaviside, είναι η κατανομή $\delta$ του Dirac.

dudx = diff(unit_step(x),x) ; print dudx

dirac_delta(x)

print type(dirac_delta)

<class 'sage.functions.generalized.FunctionDiracDelta'>

Η "συνάρτηση" $\delta$ του Dirac, στο Sage ανήκει κι αυτή στην κλάση sage.functions.generalized.FunctionUnitStep. Η βοήθεια που μας προσφέρει το Sage για την dirac_delta είναι

print sage.functions.generalized.FunctionDiracDelta.__doc__

    The Dirac delta (generalized) function, `\delta(x)` (``dirac_delta(x)``).

    INPUT:

    -  ``x`` - a real number or a symbolic expression

    DEFINITION:

    Dirac delta function `\delta(x)`, is defined in Sage as:

        `\delta(x) = 0` for real `x \ne 0` and
        `\int_{-\infty}^{\infty} \delta(x) dx = 1`

    Its alternate definition with respect to an arbitrary test
    function `f(x)` is

        `\int_{-\infty}^{\infty} f(x) \delta(x-a) dx = f(a)`

    EXAMPLES::

        sage: dirac_delta(1)
        0
        sage: dirac_delta(0)
        dirac_delta(0)
        sage: dirac_delta(x)
        dirac_delta(x)
        sage: integrate(dirac_delta(x), x, -1, 1, algorithm='sympy')
        1

    REFERENCES:

    -  http://en.wikipedia.org/wiki/Dirac_delta_function

    

Από μαθηματική σκοπιά, και οι δυο "ορισμοί" της "συνάρτησης" $\delta$ του Dirac, στο Sage, είναι τουλάχιστον τραγικοί. Ο πρώτος λέει ότι η dirac_delta(x) είναι παντού μηδέν στο $\mathbb{R}$ εκτός από το σημείο $x=0$, και έχει ολοκλήρωμα στο $\mathbb{R}$ ίσο με ένα! Πράγμα αδύνατον, αφού σύμφωνα με τα βασικά στοιχεία ολοκλήρωσης κατά Riemann (ή ακόμα και κατά Lebesgue), δυο συναρτήσεις που είναι ταυτόσημες παντού εκτός από ένα πεπερασμένο πλήθος σημείων, έχουν ακριβώς το ίδιο ολοκλήρωμα! Οπότε το ολοκλήρωμά της θα έπρεπε να είναι μηδέν κι όχι 1.
Ο δεύτερος εναλλακτικός ορισμός πάει λίγο πιο βαθιά και μιλά για ελεγκτικές συναρτήσεις (test functions). Όμως, ο προγραμματιστής του Sage που υλοποίησε την διαβολική "συνάρτηση" $\delta$, θα όφειλε να γνωρίζει ότι ούτε αυτό ισχύει, ακόμα με την έννοια των γενικευμένων συναρτήσεων, αφού η "συνάρτηση" $\delta$ του Dirac δεν μπορεί να παραχθεί από καμία τοπικά ολοκληρώσιμη συνάρτηση. Η αναφορά ιδίως στην wikipedia με αφήνει άφωνο!

Από την εμπειρία μου σε διάφορα ΣΣΥ που έχω χρησιμοποιήσει κατά καιρούς, η υλοποίηση των γενικευμένων συναρτήσεων είναι αρκετά προβληματική . Δεν αποκλείεται όμως στο άμεσο μέλλον να κατασκευασθούν αλγόριθμοι που να τις χειρίζονται σωστά από μαθηματική άποψη. Αρκεί να βρεθεί ένας προγραμματιστής που να γνωρίζει πέρα από προγραμματισμό και στοιχειώδη συναρτησιακή ανάλυση.

Για να μην τρίζουν τα κόκκαλα του Laurent Schwartz, του μαθηματικού που έθεσε το σωστό μαθηματικό πλαίσιο για τις γενικευμένες συναρτήσεις, στο μάθημα θα χρησιμοποιούμε από τις "γενικευμένες συναρτήσεις" του Sage, μόνο την αλματική συνάρτηση unit_step, για πολύ συγκεκριμένους σκοπούς. Για παράδειγμα, για να κατασκευάσουμε πιο σύνθετες ασυνεχείς συναρτήσεις.

Επιπλέον, με κατάλληλες προσαρμογές στο όρισμά της, η συνάρτηση unit_step μπορεί να λειτουργήσει σαν φίλτρο, έτσι ώστε αν μια συνθήκη $Α$ ικανοποιείται ($x\geq 1$) να πάρουμε την τιμή 1, κι αν δεν ικανοποιείται ($x< 0$), να πάρουμε την τιμή 0.

Για παράδειγμα, με την βοήθεια της unit_step μπορούμε να φτιάξουμε μια νέα συνάρτηση που έχει τρία σκαλοπάτια $$ \mathrm{f}(x) = \left\lbrace \begin{array}{cc} 1\,, & x \leq 0 \,, \\ \frac{1}{3} \,, & 0 < x \leq 1\,, \\ \frac{2}{3}\,, & 1 \leq x < 2 \,,\\ 1 \,, & x>2 \,. \end{array} \right.$$

f(x) = unit_step(x)/3 + unit_step(x-1)/3 + unit_step(x-2)/3
pp = plot(f(x),(x,-1,3),exclude=[0,1,2],thickness=2); pp.show(figsize=3)

Επιπλέον, χρησιμοποιώντας την διαφορά unit_step(x-a) - unit_step(x-b) δυο βηματικών συναρτήσεων

stepsub = plot( unit_step(x+2) - unit_step(x-1), (x,-4,3), exclude=[-2,1],thickness=2); 
stepsub.show(figsize=3)

για διάφορες τιμές των a, b, μπορούμε να φτιάξουμε (γιατί;) την τμηματικά συνεχή συνάρτηση $$ $$ $$ \mathrm{g}(x) = \left\lbrace \begin{array}{cc} 0\,, & x < -2 \,, \\ x^2 \,, & -2 \leq x < 1\,, \\ x-1\,, & 1 \leq x < 2 \,, \\ 3-x\,, & 2 \leq x < 3 \,, \\ 0\,, & x \geq 3\,. \end{array} \right.$$

g(x) = x^2*( unit_step(x+2) - unit_step(x-1) ) + \
     (x-1)*( unit_step(x-1) - unit_step(x-2) ) + \
     (3-x)*( unit_step(x-2) - unit_step(x-3) )
ppp = plot(g(x),(x,-4,5), exclude=[-2,1],thickness=2); ppp.show(figsize=4)

3.3 Δουλεύοντας με συναρτήσεις

3.3.1 Περιεκτικές λίστες

Οι περιεκτικές λίστες στην Python είναι εκφράσεις της μορφής [ f(x) for x in L ] όπου L είναι μια λίστα και f(x) μια συνάρτηση. Το αποτέλεσμα είναι μια νέα λίστα που περιέχει όλες τις τιμές της συνάρτησης f(x) στα αντίστοιχα x που ανήκουν στην λίστα L. Στο Sage μπορούμε να χρησιμοποιήσουμε τις περιεκτικές λίστες για να κατασκευάσουμε συμβολικές εκφράσεις αυθαίρετου μεγάλου μήκους και σχήματος. Για παράδειγμα, μπορούμε να κατασκευάσουμε μια ακολουθία 20 συμβολικών μεταβλητών x01,x02,...,x20.

v = [var('x' + '%02d' % k) for k in range(1,21)]
print v; print type(x13)

[x01, x02, x03, x04, x05, x06, x07, x08, x09, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20]
<type 'sage.symbolic.expression.Expression'>

Τώρα υψώνουμε στον κύβο όλες τις συμβολικές μεταβλητές της παραπάνω λίστας

C = [x^3 for x in v]
print C 

[x01^3, x02^3, x03^3, x04^3, x05^3, x06^3, x07^3, x08^3, x09^3, x10^3, x11^3, x12^3, x13^3, x14^3, x15^3, x16^3, x17^3, x18^3, x19^3, x20^3]

Για να κατασκευάσουμε μια συμβολική έκφραση ας προσθέσουμε όλους τους όρους της τελευταίας λίστας

e = sum(C)
print e; print type(e)

x01^3 + x02^3 + x03^3 + x04^3 + x05^3 + x06^3 + x07^3 + x08^3 + x09^3 + x10^3 + x11^3 + x12^3 + x13^3 + x14^3 + x15^3 + x16^3 + x17^3 + x18^3 + x19^3 + x20^3
<type 'sage.symbolic.expression.Expression'>

Για να επιστρέψουμε στην αναπαράσταση της τελευταίας συμβολικής έκφρασης με λίστα μπορούμε να χρησιμοποιήσουμε τη εντολή operands

ope = e.operands()
print ope; print type(ope)

[x01^3, x02^3, x03^3, x04^3, x05^3, x06^3, x07^3, x08^3, x09^3, x10^3, x11^3, x12^3, x13^3, x14^3, x15^3, x16^3, x17^3, x18^3, x19^3, x20^3]
<type 'list'>

Ας υποθέσουμε ότι θέλουμε να υψώσουμε στην k-οστή δύναμη την k μεταβλητή στην λίστα με το όνομα ope. Αυτό επιτυγχάνεται ως εξής

r = [ope[k]^k for k in range(len(ope))]
print r

[1, x02^3, x03^6, x04^9, x05^12, x06^15, x07^18, x08^21, x09^24, x10^27, x11^30, x12^33, x13^36, x14^39, x15^42, x16^45, x17^48, x18^51, x19^54, x20^57]

Παρατηρούμε ότι ο καθένας από τους παραπάνω όρους είναι συμβολική έκφραση

print r[15] ; print type(r[15])
print r[15].degree(x16)

x16^45
<type 'sage.symbolic.expression.Expression'>
45

οπότε με την κατάλληλη επιλογή της μεταβλητής στο όρισμα της εντολής degree() λαμβάνουμε την δύναμη 45 για την μεταβλητή x16. Με παρόμοιο τρόπο μπορούμε να πάρουμε τις δυνάμεις όλων των μεταβλητών στην λίστα

print [(r[k]).degree(v[k]) for k in range(len(ope))]

[0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51, 54, 57]

Ας υποθέσουμε ότι θέλουμε να επιλέξουμε από την λίστα εκείνους τους όρους οι οποίοι έχουν δύναμη μικρότερη από 13. Αυτό μπορεί να υπολοιηθεί με την αλματική συνάρτηση unit_step ως εξής:

print r[2], ':', unit_step(13 - r[2].degree(v[2]))
print r[8], ':', unit_step(13 - r[8].degree(v[8]))

x03^6 : 1
x09^24 : 0

Το αποτέλεσμα μας πληροφορεί ότι η μεταβλητή x03 έχει δύναμη μικρότερη από 13, αφού η συνάρτηση unit_step (με το κατάλληλο όρισμα) μας επέστρεψε 1, ενώ η x09 δεν έχει, αφού η unit_step μας έδωσε 0. Με αυτόν τον τρόπο μπορούμε να περάσουμε από ένα φίλτρο όλους τους όρους με δύναμη μεγαλύτερη από 13.

Ορίζουμε με τον τελεστή lambda της Python μια συνάρτηση που την ονομάζουμε filter και η οποία έχει δυο ορίσματα, το x και το k, και μας δίνει ως αποτέλεσμα το x πολλαπλασιασμένο με την τιμή 0 ή 1, ανάλογα αν η δύναμη της μεταβλητής x είναι μεγαλύτερη του 13 ή όχι, αντίστοιχα. Έπειτα, εφαρμόζουμε την συνάρτηση filter στην λίστα που μας ενδιαφέρει, όπου ο δείκτης k διατρέχει τα στοιχεία της λίστας r[k], και τέλος καταχωρούμε το αποτέλεσμα σε μια νέα λίστα s.

filter = lambda x, k: x*unit_step(13 - x.degree(v[k]))
s = [filter(r[k],k) for k in range(len(r))]
print s

[1, x02^3, x03^6, x04^9, x05^12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Στην νέα λίστα μπορούμε, για παράδειγμα, να αθροίσουμε τα στοιχεία της

t = sum(s)
print t

x05^12 + x04^9 + x03^6 + x02^3 + 1

ή να πάρουμε το γινόμενο των πέντε πρώτων μη μηδενικών όρων της

gin = prod(s[k] for k in (0 .. 4)) ; print gin

x02^3*x03^6*x04^9*x05^12

3.3.2 Σύνθεση συναρτήσεων

Θεωρούμε την παραβολική τροχιά ενός βλήματος που εκτοξεύθηκε από το έδαφος για $t=0\,\rm sec$ και προσέκρουσε στο έδαφος για $t=45\,\rm sec$, 12 χιλιόμετρα, πιο μακρυά. Θεωρούμε ότι το βλήμα ταξιδεύει με σταθερή οριζόντια ταχύτητα, και ονομάζουμε $f(t)$ την συνάρτηση που μας δίνει το ύψος του βλήματος σε κάθε χρονική στιγμή. Πρώτα απ' όλα, υλοποιούμε το παραβολικό σχήμα της τροχιάς του βλήματος από μια συνάρτηση $y(x)$

y(x) = x*(12 - x)
p = plot(y(x), (x, 0, 12)); p.show(figsize=3.5)

Τώρα θέλουμε την συνάρτηση $f(t)$ ως συνάρτηση του χρόνου $t$ (σε sec). Επειδή υποθέσαμε ότι το βλήμα ταξιδεύει με σταθερή οριζόντια ταχύτητα, η οριζόντια απόσταση $x$ (σε χλμ) που διανύει το βλήμα είναι μια γραμμική σχέση του $t$, και ίση με $x = \frac{12}{45}\, t$. Οπότε για $t=0$, $x=0$ και για $t=45$, $x=12$.

x(t) = 12/45*t
print 'x at 0 :', x(0)
print 'x at 45 :', x(45)

x at 0 : 0
x at 45 : 12

Συνεπώς, το ύψος $y$ του βλήματος, ως συνάρτηση του χρόνου t, είναι η σύνθεση της συνάρτησης $y(x)$ της παραβολής, με την γραμμική συνάρτηση $x=\frac{12}{45}\,t$, δηλαδή $f(t)=y(x(t))$

f(t) = y(x(t))
print 'altitude at 0 sec:', f(0)
print 'altitude at 22.5 sec:', f(22.5)
print 'altitude at 45 sec:', f(45)

altitude at 0 sec: 0
altitude at 22.5 sec: 36.0000000000000
altitude at 45 sec: 0

Με την σύνθεση συναρτήσεων έχουμε διαχωρίσει το γεωμετρικό σχήμα της τροχιάς του βλήματος από την εξέλιξη της τροχιάς στον χρόνο. Για παράδειγμα, ας μειώσουμε στο μισό την ταχύτητα του βλήματος

h(t) = t/2
hf(t) = f(h(t))
print 'altitude at 45 sec:', hf(45)
print 'altitude at 90 sec:', hf(90)

altitude at 45 sec: 36
altitude at 90 sec: 0

Παρατηρούμε ότι το βλήμα φτάνει στο μέγιστο ύψος, για $t=45\,\rm sec$, και προσκρούει στο έδαφος για $t=90\,\rm sec$. Μπορούμε να δούμε μια τρισδιάστατη απεικόνιση της τροχιάς του βλήματος με τις παρακάτω εντολές

p1 = parametric_plot3d([x(t), 0, f(t)], (t, 0, 45),thickness=4,color='blue',figsize=4); 
show(p1,viewer='tachyon')

3.3.2 Σύνθεση μιας συνάρτησης n-φορές με τον εαυτό της

Ας υποθέσουμε ότι θέλουμε με την μέθοδο του Newton να βρούμε τις ρίζες μιας εξίσωσης της μορφής $f(x)=0$

def newton_step(f, x):
    n(x) = x - f(x)/diff(f,x)
    return n

Για χάρη ευκολίας ας θεωρήσουμε ότι θέλουμε να βρούμε μια προσεγγιστική τιμή για την τετραγωνική ρίζα του 2, ως την λύση της εξίσωσης $f(x)=x^2-2=0$.

sqr(x) = x^2 - 2
our_sqrt = newton_step(sqr, x)
print our_sqrt

x |--> x - 1/2*(x^2 - 2)/x

Παρατηρούμε ότι λόγω της μορφής που έχει η συνάρτηση newton_step που ορίσαμε, αν θέσουμε ως αρχική προσέγγιση έναν ακέραιο και συνθέσουμε επαναληπτικά την συνάρτησή μας, θα πάρουμε ως αποτέλεσμα μια ρητή προσέγγιση της λύσης μας, δηλ. του $\sqrt{2}$.

print our_sqrt(our_sqrt(our_sqrt(2)))

577/408

Στο Sage υπάρχει ένας "τεμπέλικος" τρόπος για πάρουμε την σύνθεση μιας συνάρτησης με τον εαυτό της όσες φορές θέλουμε! Απλά πολλαπλασιάζουμε κατάλληλα χαρακτήρες (strings)

s = 'our_sqrt('*5 + '2' + ')'*5
print s; print type(s)

our_sqrt(our_sqrt(our_sqrt(our_sqrt(our_sqrt(2)))))
<type 'str'>

και υπολογίζουμε το αποτέλεσμα με αφετηρία ένα αρχικό σημείο, για παράδειγμα για $x=2$.

v = eval(s)
print v

886731088897/627013566048

Για να συγκρίνουμε την παραπάνω προσέγγιση με την πραγματική τιμή του $\sqrt{2}$ δίνουμε

print v.n(digits=30)
print sqrt(2).n(digits=30)

1.41421356237309504880168962350
1.41421356237309504880168872421

Φυσικά, με τον παραπάνω τρόπο μπορούμε να ορίσουμε μια νέα συμβολική έκφραση η οποία είναι η σύνθεση μιας οποιασδήποτε συνάρτηση n-φορές με τον εαυτό της! Για ευκολία παίρνουμε n=7 φορές σύνθεση την συνάρτηση του ημιτόνου με τον εαυτό του

sin_7 = 'sin('*7 + 'x' + ')'*7
print sin_7; print type(sin_7)

sin(sin(sin(sin(sin(sin(sin(x)))))))
<type 'str'>

sin_7x(x) = sin_7
print sin_7x(x); print type(sin_7x)

sin(sin(sin(sin(sin(sin(sin(x)))))))
<type 'sage.symbolic.expression.Expression'>

print sin_7x(pi/2).n(digits=30)

0.554016390755629630333389489401

3.4 Συμβολική και αριθμητική παραγώγιση

3.4.1 Συμβολική παραγώγιση

Στο Sage παραγωγίζουμε με την εντολή diff

print 'the derivative of sin(x) is', diff(sin(x),x)

the derivative of sin(x) is cos(x)

Για να παραγωγίσουμε διαδοχικά ως προς μια μεταβλητή υπάρχει ένα επιπλέον όρισμα στην εντολή diff

print [diff(sin(x), x, k) for k in range(5)]

[sin(x), cos(x), -sin(x), -cos(x), sin(x)]

Μπορούμε να παραγωγίζουμε συναρτήσεις πολλών μεταβλητών τις οποίες έχουμε ορίσει εμείς

f(x,y) = x^2*y + 2*x*y + x
print 'f =', f
print 'the derivative of f is', f.diff() ; print type(f.diff())

f = (x, y) |--> x^2*y + 2*x*y + x
the derivative of f is (x, y) |--> (2*x*y + 2*y + 1, x^2 + 2*x)
<class 'sage.modules.vector_callable_symbolic_dense.FreeModule_ambient_field_with_category.element_class'>

από όπου παρατηρούμε ότι το Sage μας επέστρεψε ένα διανυσματικό πεδίο. Οι συνιστώσες του δεν είναι άλλες από τις μερικές παραγώγους της $f(x,y)$ ως προς $x$ και $y$, αντίστοιχα. Με άλλα λόγια, αν δεν δώσουμε όρισμα στην εντολή diff για μια συνάρτηση πολλών μεταβλητών, το Sage μας επιστρέφει την κλίση (gradient) της $f$, η οποία συμβολίζεται με $${\rm{grad}} \, f \,\, \mbox{ή} \,\,\nabla f := \left( \frac{\partial f}{\partial x}, \frac{\partial f}{\partial x}\right)\,.$$ Πράγματι

print f.diff(x)
print f.diff(y)

(x, y) |--> 2*x*y + 2*y + 1
(x, y) |--> x^2 + 2*x

Τι θα συμβεί άραγε αν ζητήσουμε από το Sage να δράσει δυο φορές στην $f$ με το diff δίχως όρισμα;

H = f.diff().diff() ; print H ; print type(H)

[    (x, y) |--> 2*y (x, y) |--> 2*x + 2]
[(x, y) |--> 2*x + 2       (x, y) |--> 0]
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>

Το αποτέλεσμα είναι ένας πίνακας που περιέχει τις δεύτερες μερικές παραγώγους της $f(x,y)$, και είναι γνωστός με το όνομα Hessian. Τα στοιχεία του πίνακα είναι $H_{i\,j} = \frac{\partial^2 f}{\partial x_i \, \partial x_j}$, ή ισοδύναμα την μορφή πίνακα, για το συγκεκριμένο παράδειγμα, $$ $$ $$ \left( H_{ij}\right ) = \left( \begin{array}{cc} \frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \, \partial y} \\ \frac{\partial^2 f}{\partial y \, \partial x} & \frac{\partial^2 f}{\partial y^2} \end{array} \right)$$ $$ $$ Για να πάρουμε ένα στοιχείο από την πίνακα αυτό δίνουμε τα ονόματα των μεταβλητών στην εντολή diff. Για παράδειγμα, για να πάρουμε την μικτή μερική παράγωγο $\frac{\partial^2 f}{\partial x \, \partial y}$, δίνουμε

print f.diff(x).diff(y)

(x, y) |--> 2*x + 2

3.4.1 Αριθμητική παραγώγιση

Στο Sage η αριθμητική παραγώγιση είναι διαθέσιμη μέσω της scipy.misc

from scipy.misc import derivative as numdif
def f(x):
    return exp(x)
print 'exact derivative :', exp(1.0)
print '1st order approx :', numdif(f, 1.0, 1.0e-2)

exact derivative : 2.71828182845905
1st order approx : 2.71832713338

Η αριθμητική παραγώγιση χρησιμοποιεί κεντρικές διαφορές. Με την γνωστή μέθοδο της παρεμβολής μπορούμε να βελτιώσουμε τη αριθμητική ακρίβεια στην προσέγγιση

ND = [numdif(f, 1.0, 1.0e-2, order=k) for k in range(3, 12, 2)]
for a in ND: print a
print exp(1.0)

2.71832713338
2.71828182755
2.71828182846
2.71828182846
2.71828182846
2.71828182845905

όπου το k που δηλώνει την τάξη πρέπει να είναι περιττός θετικός ακέραιος.

3.4.1 Πεπλεγμένη παραγώγιση

Στο Sage μπορούμε να δηλώσουμε τις παραγώγους συμβολικά. Ας υποθέσουμε ότι η (εξαρτημένη) συμβολική μεταβλητή y, είναι συνάρτηση της (ανεξάρτητης) συμβολικής μεταβλητής x

var('x')
y = function('y')(x)
dy = y.diff(x)
print dy, type(dy)

diff(y(x), x) <type 'sage.symbolic.expression.Expression'>

Το παραπάνω αποτέλεσμα μπορούμε να το χρησιμοποιήσουμε για να εφαρμόσουμε πεπλεγμένη παραγώγιση

circle = x^2 + y^2 - 1 == 0
print circle; print type(circle)

x^2 + y(x)^2 - 1 == 0
<type 'sage.symbolic.expression.Expression'>

Η παραπάνω εξίσωση του κύκλου ορίζει, με έμμεσο τρόπο, πώς η συμβολική μεταβλητή $y$ εξαρτάται από το $x$. Παραγωγίζοντας την εξίσωση του κύκλου ως προς $x$, παίρνουμε την εξίσωση

dc = circle.diff(x)
print dc

2*y(x)*diff(y(x), x) + 2*x == 0

την οποία φυσικά μπορούμε να την λύσουμε ως προς την παράγωγο $y'(x)$, που δεν είναι άλλη από την συμβολική μεταβλητή που ονομάσαμε dy

yx = solve(dc, dy)
print yx

[
diff(y(x), x) == -x/y(x)
]

3.5 Ολοκληρώματα κι αθροίσματα

3.5.1 Ορισμένα κι αόριστα ολοκληρώματα

Το θεµελιώδες θεώρηµα του απειροστικού λογισµού μας πληροφορεί ότι η αντίστροφη διαδικασία της παραγώγισης είναι η ολοκλήρωση. Δοσμένης μιας παραγωγίσιμης συνάρτησης μπορούμε, κατ'αρχήν, να βρούμε την συνάρτηση της παραγώγου της σε κλειστή μορφή. Όμως, γενικά, δεν ισχύει το ίδιο για το ολοκλήρωμά της. Μάλιστα ένα αόριστο ολοκλήρωμα, ως συνάρτηση του άνω ορίου, αποτελεί έναν εναλλακτικό ορισμό μιας νέας πιο σύνθετης συνάρτησης. Τα κλασικά παραδείγματα είναι οι εναλλακτικοί ορισμοί των συναρτήσεων του λογαρίθμου και του τόξου εφαπτομένης: $$ \log x = \int_1^x \frac{1}{t} \, dt \,,\qquad \arctan x = \int_0^x \frac{1}{1+t^2}\,dt\,.$$ Ας θεωρήσουμε την συνάρτηση $f(x) = e^{-x^2}$, της οποίας η παράγωγος είναι

f(x) = exp(-x^2)
df = f.diff(x); print df

x |--> -2*x*e^(-x^2)

Η γραφική παράσταση της $f(x) = e^{-x^2}$ μοιάζει με μια καμπάνα γύρω από το $x=0$, η οποία σβήνει στο μηδέν πάρα πολύ γρήγορα.

fplot = plot(f(x),(x,-3,3),figsize=4) ; fplot.show()

Ας ζητήσουμε από το Sage το ορισμένο ολοκλήρωμα της $f(x) = e^{-x^2}$, στο διάστημα $(-\infty, \infty)$

orismeno = integral(f(x),(x,-oo,oo)); print orismeno

sqrt(pi)

δηλαδή $$ \int_{-\infty}^\infty e^{-t^2}\,dt = \sqrt{\pi}\,. $$ Παρατηρήστε ότι στο Sage το άπειρο ($\infty$) δηλώνεται με δυο κολλητά όμικρον (οο). Επειδή η $f(x) = e^{-x^2}$ είναι άρτια συνάρτηση σε όλο το πεδίο ορισμού της, δηλαδή το $\mathbb{R}$, ισχύει ότι $$ \int_{0}^\infty e^{-t^2}\,dt = \frac{\sqrt{\pi}}{2} \qquad \Rightarrow \qquad \frac{2}{\sqrt{\pi}} \, \int_{0}^\infty e^{-t^2}\,dt =1 \,. $$

Χρησιμοποιώντας την γνωστή ιδιότητα της παρεμβολής του ολοκληρώματος Riemann, για οποιοδήποτε $x\in \mathbb{R}$, (γιατί;) έχουμε ότι $$\frac{2}{\sqrt{\pi}} \, \left( \int_{0}^x \,\,+ \int_{x}^\infty e^{-t^2}\,dt \right) =1 \quad \Rightarrow \quad \frac{2}{\sqrt{\pi}} \, \int_{0}^x e^{-t^2}\,dt + \frac{2}{\sqrt{\pi}} \, \int_{x}^\infty e^{-t^2}\,dt = 1 \qquad (*) $$

Στην τελευταία ισότητα το πρώτο ολοκλήρωμα είναι ένα αόριστο ολοκλήρωμα της συνάρτησης $f(x) = e^{-x^2}$, πολλαπλασιασμένο με $\frac{2}{\sqrt{\pi}}$. Ας ρωτήσουμε το Sage ποιό είναι το αόριστο ολοκλήρωμα με την εντολή εντολή integral, η οποία δέχεται για ορίσματα την υπό ολοκλήρωση συνάρτηση και την μεταβλητή στην οποία ολοκληρώνουμε

ex2 = exp(-x^2)
aoristo = 2/sqrt(pi)*integral(ex2,x); print aoristo

erf(x)

Το αποτέλεσμα είναι μια συνάρτηση που το Sage την ονομάζει erf(x), και είναι γνωστή ως η συνάρτηση σφάλματος (error function). Αν πάρουμε βοήθεια για την συνάρτηση αυτή από το Sage μας πληροφορεί ότι ορίζεται ακριβώς με την μορφή του αόριστου ολοκληρώματος

$$\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt\,, $$

σε απόλυτη συμφωνία με τους προηγούμενους υπολογισμούς μας. Η γραφική παράσταση της συνάρτησης σφάλματος είναι

erfplot = plot(erf(x),(x,-7,7),figsize = 4) ; erfplot.show()

Επιπλέον, από την σχέση $(*)$ παίρνουμε ότι $$ \frac{2}{\sqrt{\pi}} \, \int_{x}^\infty e^{-t^2}\,dt = 1 -{\rm{erf}}(x) $$

Η συνάρτηση $\rm{erfc}(x) = 1-\rm{erf}(x)$, ή ισοδύναμα ${\rm{erfc}}(x) = \frac{2}{\sqrt{\pi}} \, \int_{x}^\infty e^{-t^2}\,dt$, ονομάζεται η συμπληρωματική συνάρτηση σφάλματος $\rm{erfc}(x) = 1 - \rm{erf}(x)$, και η γραφική της παράσταση είναι

erfcplot = plot(1-erf(x),(x,-7,7),figsize = 4) ; erfcplot.show()

Στο Sage δεν έχει κάθε συμβολική συνάρτηση ορισμένο ή αόριστο ολοκλήρωμα. Για παράδειγμα

r = integral(tan(exp(-x^2)),(x,0,1)) ; print r

integrate(tan(e^(-x^2)), x, 0, 1)

Αν η συνάρτηση της οποίας αναζητούμε το ορισμένο ολοκλήρωμα σε κάποιο διάστημα είναι ολοκληρώσιμη κατά Riemann, μπορούμε να ζητήσουμε από το Sage μια αριθμητική προσέγγισή του

print r.n()

0.996025062538

3.5.2 Βοηθώντας το Sage

Ας ζητήσουμε από το Sage να μας δώσει το ολοκλήρωμα Riemann της $f(x)=x^2$, όπου τα άκρα της ολοκλήρωσης είναι συμβολικές εκφράσεις

a , b = var(' a , b')
print integral(x^2 , (x,a,b))

-1/3*a^3 + 1/3*b^3

Το αποτέλεσμα είναι αναμενόμενο αφού γνωρίζουμε ότι το ολοκλήρωμα είναι $$ $$ $$ \int_a^b x^2\, dx = \frac{1}{3}\, \left. x^3\right|_a^b = \frac{1}{3}\,b^3-\frac{1}{3}\,a^3 \,,$$ $$ $$ το οποίο γεωμετρικά παριστάνει το εμβαδό ανάμεσα στον x-άξονα, την παραβολή $y=x^2$, και τις κατακόρυφες ευθείες $x=a$, $x=b$, όπως δείχνεται στο παρακάτω σχήμα

f(x) = x^2
fplot = plot(f(x),(x,-2,3), ticks=[[-2.0,3.0],[4,9]], \
             tick_formatter=[["$a$","$b$"],["$f(a)$","$f(b)$"]],\
             fill='axis',figsize=3); fplot.show()

Ας ζητήσουμε τώρα από το Sage το ολοκλήρωμα της $y=\frac{1}{x^2}$ στο διάστημα $(a,b)$

a , b = var(' a , b')
pprint integral(1/x^2 , (x,a,b))

  File "<ipython-input-62-958958d8eff2>", line 2
       
ValueError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(b>0)', see 'assume?' for more details)
Is b positive, negative or zero? 

Το Sage μας επιστρέφει λάθος (ακριβέστερα το Maxima το οποίο εκτελεί την συμβολική ολοκλήρωση), κι ύστερα από αρκετές γραμμές κώδικα κατέληξε στο συμπέρασμα που σημειώνει παραπάνω με κόκκινο χρώμα

Με μια πιο προσεχτική ματιά, η συνάρτηση της οποίας ζητάμε το ολοκλήρωμα Riemann δεν ορίζεται στο σημείο $x=0$, και καθώς το $x$ πλησιάζει το μηδέν από αριστερά ή από δεξιά απειρίζεται θετικά, δηλαδή η συνάρτησή μας έχει ασύμπτωτο τον $y$ άξονα. Έτσι το Sage έχει κάθε λόγο να παραπονιέται, γιατί η $y=\frac{1}{x^2}$ είναι συνεχής μόνο κατά διαστήματα, στα $(-\infty,0)$ και $(0,\infty)$, και το θεμελιώδες θεώρημα του απειροστικού λογισμού δεν ισχύει σε ένωση ξένων διαστημάτων. Η γραφική παράσταση της συνάρτησής μας είναι

f(x) = 1/x^2
fplot = plot(f(x),(x,-2,2),ymin=1,ymax=50,figsize=4,detect_poles=True)
fplot.set_legend_options(back_color=(1,0.9,0.9), shadow=False); fplot.show()

Οπότε για να πάρουμε ένα συνεπές αποτέλεσμα πρέπει να περιορίσουμε κατάλληλα τις τιμές των συμβολικών μεταβλητών $a,b$, στα άκρα της ολοκλήρωσης. Για παράδειγμα, $b>a>0$.

assume( a > 0, b > a)
print integral(1/x^2, (x,a,b))

1/a - 1/b

Σε κάθε στιγμή μπορούμε να πληροφορηθούμε για τις υποθέσεις που έχουμε επιβάλλει στις συμβολικές μεταβλητές μας με την εντολή assumptions

print assumptions()

[a > 0, b > a]

Οι περιορισμοί που επιβάλλουμε στις συμβολικές μεταβλητές την εντολή assumptions αναιρούνται με την εντολή forget

forget(b > a); print assumptions()
forget(a > 0); print assumptions()

[a > 0]
[]

3.5.2 Συμβολικά αθροίσματα

Το Sage μπορεί να βρει σε κλειστή μορφή το άθροισμα συμβολικών εκφράσεων. Για παράδειγμα ας θεωρήσουμε το άθροισμα των n πρώτων όρων της αριθμητικής προόδου $$ \sum_{k=1}^n k = 1+2+3+\cdots + n\,.$$

var('k, n')
s = sum(k, k, 1, n)
s.show()
f = s.factor()
f.show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, n^{2} + \frac{1}{2} \, n$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, {\left(n + 1\right)} n$$

Το ενδιαφέρον είναι ότι αν μια σειρά $\sum_{k=1}^\infty a_k$ συγκλίνει, το Sage είναι ικανό να υπολογίσει το άθροισμα της σειρά. Για παράδειγμα, ας θεωρήσουμε την σειρά $ \sum_{k=1}^\infty \frac{1}{k^2}$, η οποία από τα κριτήρια σύγκλισης σειρών που μάθαμε στον απειροστικό λογισμό, γνωρίζουμε ότι συγκλίνει. Σε ποιόν αριθμό να συγκλίνει άραγε;

sum(1/k^2, k, 1, oo).show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{6} \, \pi^{2}$$

Ας πάρουμε τώρα, την εναλλασσόμενη σειρά $$ \sum_{k=1}^\infty (-1)^k \frac{1}{k}\,,$$ η οποία, από τα κριτήρια των εναλλασσόμενων σειρών, γνωρίζουμε ότι και αυτή συγκλίνει.

print sum((-1)^(k)/k, k, 1, oo)

-log(2)

Χρησιμοποιώντας τις περιεκτικές λίστες της Python μπορούμε να ελέγξουμε το τελευταίο αποτέλεσμα που μας δίνει το Sage

L = [float((-1)^k/k) for k in range(1, 100001)]
print -sum(L)
print (log(2)).n(digits=16)

0.693142180585
0.6931471805599453

απ' όπου συμπεραίνουμε ότι η σειρά συγκλίνει πολύ αργά.

Για τον υπολογισμό του αθροίσματος της σειράς $\sum_{n=0}^\infty r^n$, το Sage χρειάζεται την βοήθειά μας για να περιορίσουμε τις τιμές της συμβολικής έκφρασης $r$, σε $|r|<1$,

var('r , n')
assume( abs(r) < 1)
print sum(r^n, n,0,oo)

-1/(r - 1)

3.6 Σειρές Taylor, δυναμοσειρές, προσέγγιση και όρια

3.6.1 Σειρές Taylor

Στο Sage το ανάπτυγμα μιας συνάρτησης γύρω από ένα σημείο δίνεται με την εντολή taylor. Ας βρούμε την σειρά Taylor της συνάρτησης cos(x) γύρω από το x=0, μέχρι όρους τάξης 6

Tcos = taylor(cos(x),x,0,6); print Tcos; print type(Tcos)

-1/720*x^6 + 1/24*x^4 - 1/2*x^2 + 1
<type 'sage.symbolic.expression.Expression'>

Το αποτέλεσμα δεν είναι μια σειρά αλλά το πολυώνυμο που προκύπτει αν αποκόψουμε από την σειρά όλους τους όρους με τάξη μεγαλύτερη από 6. Επίσης παρατηρούμε ότι το αποτέλεσμα είναι μια συμβολική έκφραση. Μπορούμε να ζητήσουμε να μας δώσει το αποτέλεσμα καθαρά συμβολικά

var('x, a')
f = function('f')(x)
Tf = taylor(f(x), x, a, 4)
print Tf; print type(Tf)

1/24*(a - x)^4*diff(f(a), a, a, a, a) - 1/6*(a - x)^3*diff(f(a), a, a, a) + 1/2*(a - x)^2*diff(f(a), a, a) - (a - x)*diff(f(a), a) + f(a)
<type 'sage.symbolic.expression.Expression'>

Σαν ένα άλλο παράδειγμα θεωρούμε την ακολοθουθία των αριθμών Fibonacci η οποία εμφανίζεται ως οι συντελεστές της σειράς Taylor γύρω από το x=0, της ρητής συνάρτησης $$ g(x) = \frac{1}{1-x-x^2}\,.$$ Ας δούμε τους 10 πρώτους αριθμούς της ακολουθίας αριθμών Fibonacci

g(x)  = 1/(1-x-x^2)
Tg = taylor(g(x), x, 0, 20)
print Tg;
print Tg.coefficients()

10946*x^20 + 6765*x^19 + 4181*x^18 + 2584*x^17 + 1597*x^16 + 987*x^15 + 610*x^14 + 377*x^13 + 233*x^12 + 144*x^11 + 89*x^10 + 55*x^9 + 34*x^8 + 21*x^7 + 13*x^6 + 8*x^5 + 5*x^4 + 3*x^3 + 2*x^2 + x + 1
[[1, 0], [1, 1], [2, 2], [3, 3], [5, 4], [8, 5], [13, 6], [21, 7], [34, 8], [55, 9], [89, 10], [144, 11], [233, 12], [377, 13], [610, 14], [987, 15], [1597, 16], [2584, 17], [4181, 18], [6765, 19], [10946, 20]]

Το Sage μπορεί να μας δώσει πολυωνυμικά αναπτύγματα σε σειρά Taylor και για συναρτήσεις πολλών μεταβλητών

var('x,y')
h = cos(x)*sin(y)
print taylor(h,(x,0),(y,pi/2),5)

1/384*(pi - 2*y)^4 + 1/16*(pi - 2*y)^2*x^2 + 1/24*x^4 - 1/8*(pi - 2*y)^2 - 1/2*x^2 + 1

3.6.2 Σειρές Taylor και SymPy

Για να δούμε ποιά συστήματα χρησιμοποιεί το Sage για τον υπολογισμό των σειρών Taylor δίνουμε

from sage.misc.citation import get_systems
print get_systems('taylor(cos(x), x, 0, 6)')

['ginac', 'Maxima']

Μπορούμε εναλλακτικά να χρησιμοποιήσουμε το SymPy που μας δίνει το αποτέλεσμα στην μορφή Python

from sympy.series import series
print series(cos(x), x0=0, n=10)

1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)

Αν δώσουμε n=None, ως αποτέλεσμα θα πάρουμε μια γεννήτρια συμβολικών εκφράσεων με την βοήθεια της οποίας μπορούμε να βρούμε επαναληπτικά τον επόμενο όρο της σειράς Taylor με την εντολή next

g = series(cos(x), x0=0, n=None)
print g; print type(g)

<generator object <genexpr> at 0x7f230d8cec80>
<type 'generator'>

print g.next()
print g.next()

1
-x**2/2

Ακόμα καλύτερα, χρησιμοποιώντας και πάλι την περικτική λίστα της Python [ F for k in L ] παίρνουμε μια λίστα σε μορφή Python, όπου τα στοιχεία της είναι οι όροι της σειράς Taylor,

print [g.next() for k in range(5)]

[x**4/24, -x**6/720, x**8/40320, -x**10/3628800, x**12/479001600]

3.6.3 Δυναμοσειρές

Όπως είδαμε στο Sage οι "σειρές" Taylor δεν είναι δυναμοσειρές γύρω από κάποιο σημείο $x_0$, αλλά πολυώνυμα που προκύπτουν αν αποκόψουμε από την σειρά όλους τους όρους της από μια τάξη και πάνω, ανάλογα με την αντίστοιχη τάξη της παραγώγου. Όμως το Sage μας παρέχει την δυνατότητα να μετατρέψουμε ένα πολυώνυμο (γύρω από το μηδέν) σε μια δυναμοσειρά

var('x')
q  = x^3+5*x + 7
pq = q.power_series(QQ)
print pq; print type(pq)

7 + 5*x + x^3 + O(x^4)
<type 'sage.rings.power_series_poly.PowerSeries_poly'>

Ο λόγος που ίσως μας είναι χρήσιμη η μετατροπή ενός πολυωνύμου σε μια δυναμοσειρά, είναι ότι για τις δυναμοσειρές μπορεί να ορισθεί η διαίρεση (μέχρι όρους τάξης n) και συνεπώς μια δυναμοσειρά έχει έναν αντίστροφο (μέχρι όρους τάξης n), το οποίο δεν ισχύει για πολυώνυμα. Ας δούμε ένα παράδειγμα

invpq = 1/pq
print invpq; print type(invpq); print pq*invpq

1/7 - 5/49*x + 25/343*x^2 - 174/2401*x^3 + O(x^4)
<type 'sage.rings.power_series_poly.PowerSeries_poly'>
1 + O(x^4)

Παρατηρούμε ότι το γινόμενο της δυναμοσειράς pq με την invpq μας δίνει 1 μέχρι όρους τάξης 4. Η τάξη αυτή είναι η τάξη της ακρίβειας των υπολογισμών και δεν είναι ίση με τον βαθμό του πολυωνύμου που εμφανίζεται στην δυναμοσειρά

print pq.prec(); print invpq.prec(); 

4
4

print pq.degree()

3

Μπορούμε επίσης να μετρατρέψουμε μια δυναμοσειράς σε μια δυναμοσειρά μκρότερης ακρίβειας χρησιμοποιώντας την μέθοδο trancate_powerseries()

trinvpq = invpq.truncate_powerseries(2); 
print trinvpq; print type(trinvpq)

1/7 - 5/49*x + O(x^2)
<type 'sage.rings.power_series_poly.PowerSeries_poly'>

Μπορούμε επιπλέον να αποκόψουμε μια δυναμοσειρά σε ένα πολυώνυμο με την μέθοδο polynomial(), ή εναλλακτικά με την μέθοδο trancate()

print invpq.polynomial()
print invpq.truncate(2)

-174/2401*x^3 + 25/343*x^2 - 5/49*x + 1/7
-5/49*x + 1/7

3.6.3 Προσεγγίσεις

Με την προσέγγιση Padé προσεγγίζουμε μια συνάρτηση από μια ρητή συνάρτηση, όπου ο βαθμός του πολυωνύμου του αριθμητή και του παρανομαστή είναι δοσμένοι θετικοί ακέραιοι. Οι συντελεστές των πολυωνύμων είναι ρητοί αριθμοί και προσδιορίζονται από την μέθοδο Padé, στην οποία δεν θα επεκταθούμε περισσότερο εδώ. Συνήθως η προσέγγιση Padé μας δίνει μια καλύτερη προσέγγιση από τις σειρές Taylor, ιδίως όταν η συνάρτηση που θέλουμε να προσεγγίσουμε έχει πόλους. Στο Sage η προσέγγιση Padé μιας συνάρτησης υπολοιείται με την "μέθοδο" pade, που δέχεται σαν ορίσματα τον βαθμό του αριθμητή και του παρανομαστή.

z = PowerSeriesRing(QQ, 'z').gen()
pad = cos(z).pade(5,6)
print pad; print type(pad)

(60514/59*z^4 - 1428000/59*z^2 + 3155040/59)/(z^6 + 3814/59*z^4 + 149520/59*z^2 + 3155040/59)
<class 'sage.rings.fraction_field_element.FractionFieldElement_1poly_field'>

Για να δούμε πόσο καλή είναι η προσέγγιση Padé της συνάρτησης του συνημιτόνου ας υπολογίσουμε την τιμή της σε ένα σημείο κι ας την συγκρίνουμε με την πραγματική τιμή στο ίδιο σημείο.

print pad(pi).n(digits=10)
print cos(pi)

-0.9970661822
-1

pp = plot(pad,(z,-2*pi,2*pi),legend_label='pad $\cos \,x$',color='red');
pc = plot(cos(x),(x,-2*pi,2*pi),legend_label='$\cos\, x$');
(pp+pc).show(figsize=4)

Θα έχουμε την ευκαιρία κατά την διάρκεια των μαθημάτων να δούμε προσεγγίσεις μιας συνάρτησης μέσω άλλων συναρτήσεων, όπως τριγωνομετρικές, όταν θα αναφερθούμε στις σειρές Fourier, και τα πολυώνυμα Chebyshev.

3.6.3 Όρια

Στο Sage μπορούμε να βρούμε όρια συναρτήσεων (συμπεριλαμβανομένων πλευρικών ορίων) με την εντολή limit, ή ισοδύναμα με την συντομογραφία της lim

var('a,n,x')
print lim( (1+a/n)^n,n=oo)

e^a

print lim(1/x - log(1/x) , x=oo ,dir='+')

+Infinity

print lim( x^(1/x) , x=oo)

1

print lim( 1/x , x=0 , dir='-')

-Infinity

print lim(sin(1/x), x = 0)

ind

Με ind το Sage εννοεί απροσδιοριστία, δηλαδή δεν υπάρχει όριο αλλά οι τιμές της συνάρτησης παραμένουν φραγμένες.

3.6 Δεσμευμένα ακρότατα

Στην τελευταία αυτή παράγραφο του κεφαλαίου θα δούμε ένα παράδειγμα που είναι ενδεικτικό της αποτελεσματικής εφαρμογής των συμβολικών και αριθμητικών εργαλείων του Sage. Υποθέτω ότι είστε εξοικειωμένοι με το πρόβλημα βελτιστοποίησης των τιμών μιας συνάρτησης πολλών μεταβλητών υπό δέσμευση, με την μέθοδο των πολλαπλασιαστών Lagrange. Πιο συγκεκριμένα.
Μας ενδιαφέρει να βρούμε τις βέλτιστες (μέγιστες κι ελάχιστες) τιμές της συνάρτησης $$ $$ $$ f(x,y,x) = x \,y\, z \,,$$ $$ $$ δοσμένου ότι τα σημεία $(x,y,z)$ είναι σημεία σφαίρας ακτίνας μονάδας και κέντρου $(0,0,0)$ $$ $$ $$ g(x,y,z)=x^2+y^2+z^2=1\,.$$ $$ $$ Στο Sage ορίζουμε με συμβολικές εκφράσεις την συνάρτηση που θέλουμε να βελτιστοποιήσουμε και την εξίσωση της δέσμευσης

var('x,y,z')
target = x*y*z
constraint = x^2 + y^2 + z^2 - 1 == 0
print 'optimize', target, 'subject to', constraint

optimize x*y*z subject to x^2 + y^2 + z^2 - 1 == 0

Για να εφαρμόσουμε την μέθοδο των πολλαπλασιαστών Lagrange, πρέπει να βρούμε τις κλίσεις των συναρτήσεων $f\,,g\,,$ οπότε πρέπει να μετατρέψουμε τις συμβολικές εκφράσεις μας σε συναρτήσεις για να είναι σε θέση το Sage να υπολογίσει συμβολικά μερικές παραγώγους

f(x,y,z) = target
g(x,y,z) = constraint.lhs()
print 'target as function :', f
print 'constraint as function :', g
gradf = f.diff()
gradg = g.diff()
print 'gradient of target :', gradf
print 'gradient of constraint :', gradg

target as function : (x, y, z) |--> x*y*z
constraint as function : (x, y, z) |--> x^2 + y^2 + z^2 - 1
gradient of target : (x, y, z) |--> (y*z, x*z, x*y)
gradient of constraint : (x, y, z) |--> (2*x, 2*y, 2*z)

Αφού υπολογίσαμε τις κλίσεις είμαστε έτοιμοι να φτιάξουμε το σύστημα των εξισώσεων που η λύση τους θα μας δώσει τα δεσμευμένα ακρότατα. Συνήθως στην μέθοδο με τους πολλαπλασιαστές Lagrange η πραγματική παράμετρος δηλώνεται με $\lambda$. Στο παρόν παράδειγμα την δηλώνουμε με την συμβολική μεταβλητή ell

var('ell')
eqs = [gradf(x,y,z)[k] - ell*gradg(x,y,z)[k] == 0 for k in range(3)]
eqs.append(constraint)
print eqs

[-2*ell*x + y*z == 0, -2*ell*y + x*z == 0, x*y - 2*ell*z == 0, x^2 + y^2 + z^2 - 1 == 0]

Τώρα δεν έχουμε παρά να λύσουμε το σύστημα των εξισώσεων για να βρούμε τα δεσμευμένα ακρότατα

sols = solve(eqs,[x,y,z,ell])
print sols
print 'found', len(sols), 'solutions'

[
[x == 0, y == -1, z == 0, ell == 0],
[x == 0, y == 1, z == 0, ell == 0],
[x == 0, y == 0, z == -1, ell == 0],
[x == 0, y == 0, z == 1, ell == 0],
[x == 1, y == 0, z == 0, ell == 0],
[x == -1, y == 0, z == 0, ell == 0],
[x == -1/3*sqrt(3), y == 1/3*sqrt(3), z == -1/3*sqrt(3), ell == 1/6*sqrt(3)],
[x == 1/3*sqrt(3), y == -1/3*sqrt(3), z == -1/3*sqrt(3), ell == 1/6*sqrt(3)],
[x == -1/3*sqrt(3), y == -1/3*sqrt(3), z == 1/3*sqrt(3), ell == 1/6*sqrt(3)],
[x == 1/3*sqrt(3), y == 1/3*sqrt(3), z == 1/3*sqrt(3), ell == 1/6*sqrt(3)],
[x == -1/3*sqrt(3), y == 1/3*sqrt(3), z == 1/3*sqrt(3), ell == -1/6*sqrt(3)],
[x == 1/3*sqrt(3), y == -1/3*sqrt(3), z == 1/3*sqrt(3), ell == -1/6*sqrt(3)],
[x == -1/3*sqrt(3), y == -1/3*sqrt(3), z == -1/3*sqrt(3), ell == -1/6*sqrt(3)],
[x == 1/3*sqrt(3), y == 1/3*sqrt(3), z == -1/3*sqrt(3), ell == -1/6*sqrt(3)]
]
found 14 solutions

Το Sage βρήκε 14 λύσεις, οπότε 14 πιθανά δεσμευμένα ακρότατα. Οι αντίστοιχες τιμές της συνάρτησης $f(x,y,z)$ που θέλουμε να βελτιστοποιήσουμε είναι

print [target.substitute(sols[k]) for k in range(0,len(sols))]

[0, 0, 0, 0, 0, 0, 1/9*sqrt(3), 1/9*sqrt(3), 1/9*sqrt(3), 1/9*sqrt(3), -1/9*sqrt(3), -1/9*sqrt(3), -1/9*sqrt(3), -1/9*sqrt(3)]

οπότε τα σημεία

$\left(\frac{-\sqrt{3}}{3}\,, \frac{\sqrt{3}}{3} ,\,\frac{-\sqrt{3}}{3}\right)$, $\left(\frac{\sqrt{3}}{3}\,, \frac{-\sqrt{3}}{3} ,\,\frac{-\sqrt{3}}{3}\right)$, $\left(\frac{\sqrt{3}}{3}\,, \frac{\sqrt{3}}{3} ,\,\frac{\sqrt{3}}{3}\right)$, $\left(\frac{-\sqrt{3}}{3}\,, \frac{-\sqrt{3}}{3} ,\,\frac{\sqrt{3}}{3}\right)$ είναι δεσμευμένα μέγιστα,

και τα σημεία

$\left(\frac{-\sqrt{3}}{3}\,, \frac{-\sqrt{3}}{3} ,\,\frac{-\sqrt{3}}{3}\right)$, $\left(\frac{\sqrt{3}}{3}\,, \frac{\sqrt{3}}{3} ,\,\frac{-\sqrt{3}}{3}\right)$, $\left(\frac{-\sqrt{3}}{3}\,, \frac{\sqrt{3}}{3} ,\,\frac{\sqrt{3}}{3}\right)$, $\left(\frac{\sqrt{3}}{3}\,, \frac{-\sqrt{3}}{3} ,\,\frac{\sqrt{3}}{3}\right)$ είναι δεσμευμένα ελάχιστα.