6. Συμβολική κι αριθμητική επίλυση εξισώσεων κι ανισώσεων¶
Πολλά προβλήματα των μαθηματικών ανάγονται στην επίλυση ενός συστήματος εξισώσεων οι οποίες μπορεί να είναι αλγεβρικές, υπερβατικές, διαφορικές, ολοκληρωτικές κτλ. Σε αυτό το κεφάλαιο θα μελετήσουμε την επίλυση συστημάτων μη-γραμμικών αλγεβρικών και υπερβατικών εξισώσεων με την βοήθεια του Sage, τόσο συμβολικά όσο και αριθμητικά.
6.1 Συμβολική επίλυση εξισώσεων¶
Ας θεωρήσουμε μια απλή τριγωνομετρική εξίσωση της μορφής $$ \sin x= \cos x$$ και ας ζητήσουμε από το Sage να μας λύσει την εξίσωση ως προς $x$
eq = sin(x)-cos(x)==0
print solve(eq,x)
Το αποτέλεσμα δεν είναι και τόσο διαφωτιστικό αφού το Sage μας επέστρεψε την εξίσωση που θέλουμε να επιλύσουμε. Ας δούμε πιo προσεκτικά ποιά συστήματα χρησιμοποιεί το Sage για την επίλυση εξισώσεων
from sage.misc.citation import get_systems
print get_systems('solve(x^2-x==0,x)')
Συνεπώς για επίλυση εξισώσεων το Sage χρησιμοποιεί το Maxima. Από την βοήθεια για την εντολή solve παρατηρούμε ότι υπάρχουν διαθέσιμες κάποιες επιλογές όταν θέλουμε να λύσουμε εξισώσεις με πολλές ρίζες, όπως οι τριγωνομετρικές εξισώσεις. Η επιλογή που μας προτρέπει να χρησιμοποιήσουμε το Sage είναι to_poly_solve = 'force'
print solve(eq,x,to_poly_solve='force')
Το αποτέλεσμα εμφανίζει δυο ελεύθερες μεταβλητές οι οποίες προέκυψαν από την διαδικασία επίλυσης που εκτελεί το Maxima. Το όνομά τους αρχίζει με το γράμμα z που σημαίνει ότι πρόκειται για μεταβλητές που διατρέχουν το σύνολο $\mathbb{Z}$ των ακεραίων. Οπότε οι λύσεις που μας δίνει το Sage είναι
$$x=-\frac{3\,\pi}{4} + 2\,k\,\pi \,, \qquad x=\frac{\pi}{4} + 2\,k\,\pi \,,\qquad k\in \mathbb{Z}\,,$$
όπου τα πολλαπλάσια του $\pi$ είναι άρτια. Όμως παρατηρούμε ότι μπορούμε να γράψουμε και τις δυο λύσεις σε μία γιατί
$$x=-\frac{3\,\pi}{4} + 2\,k\,\pi = \frac{\pi}{4} - \pi +2\,k\,\pi = \frac{\pi}{4} +(2\,k-1)\,\pi \,.$$
Οπότε, ουσιαστικά το αποτέλεσμα που μας δίνει το Sage είναι
$$ x= \frac{\pi}{4} + k\,\pi \,,\qquad k\in \mathbb{Z} \,, $$
όπου αν το $k$ είναι άρτιος έχουμε την δεύτερη λύση, ενώ αν το $k$ είναι περιττός έχουμε την πρώτη λύση. Το αποτέλεσμα μπορούμε να το ελέγξουμε άμεσα δίνοντας
var('k'); assume(k,'integer');
print eq(x = pi/4 + k*pi).trig_expand();
forget(k,'integer')
Με παρόμοιο τρόπο το Sage επιλύει συμβολικά συστήματα εξισώσεων. Για παράδειγμα
var('x y')
print solve([cos(x)*sin(y) == 0, x + y==0],x,y,to_poly_solve='force')
όπου και πάλι οι ελεύθερες μεταβλητές διατρέχουν το σύνολο $\mathbb{Z}$ των ακεραίων, γιατί το όνομά τους αρχίζει από z. Ας επιλύσουμε το ακόλουθο σύστημα εξισώσεων
var('x y')
print solve([x+y == 3, 2*x+2*y == 6],x,y)
Εδώ η ελεύθερη μεταβλητή διατρέχει το σύνολο $\mathbb{R}$ των πραγματικών αριθμών, αφού το όνομά της αρχίζει από r. Όταν είναι εφικτό, το Sage επιλύει συμβολικά, δηλαδή ακριβώς, και συστήματα εξισώσεων όπου τις εξισώσεις τις δηλώνουμε σε λίστα. Για παράδειγμα
var('x y' )
print solve([sqrt(x) + sqrt(y) == 5, x + y == 10], x, y)
6.2 Αριθμητική επίλυση εξισώσεων¶
6.2.1 Βαθμωτές υπερβατικές εξισώσεις¶
Ας θεωρήσουμε την εξίσωση
$$ \sin x + x\, \cos x =0\,,$$
και ας ζητήσουμε από το Sage να μας λύσει την εξίσωση ως προς $x$
var('x')
eq = sin(x) + x*cos(x) == 0
print solve(eq,x)
Το αποτέλεσμα δεν είναι αυτό που περιμέναμε. Το Sage διαίρεσε(!) με το $\cos x$ και μας επέστρεψε ουσιαστικά την ίδια εξίσωση της οποίας θέλουμε να βρούμε τις ρίζες. Με μια δεύτερη σκέψη, το αποτέλεσμα δεν πρέπει να μας δημιουργεί και τόσο μεγάλη εντύπωση αφού δεν αναμένουμε να υπάρχει κλειστός τύπος που να μας δίνει τις ρίζες της εξίσωσης.
Όμως, το αποτέλεσμα όπως μας το επιστρέφει το Sage μας δίνει την εξής ιδέα. Οι ρίζες της εξίσωσης είναι τα σημεία τομής του γραφήματος της συνάρτησης $y=-\tan x$ με την ευθεία $y=x$. Προτού κάνουμε μια γραφική παράσταση και των δυο συναρτήσεων ας βεβαιωθούμε ότι τα σημεία που μηδενίζεται η $\cos x$ δεν είναι ρίζες της εξίσωσής μας.
r = solve(cos(x)==0,x,to_poly_solve='force'); print r
var('k')
assume(k,'integer')
assume(k,'odd')
print eq(x=1/2*pi+k*pi).lhs().trig_expand()
forget(k,'odd'); assume(k,'even')
print eq(x=1/2*pi+k*pi).lhs().trig_expand()
forget(k,'even'); forget(k,'integer');
Πράγματι, αν $k$ είναι περιττός το αριστερό μέλος της εξίσωσης είναι -1, ενώ αν $k$ άρτιος το αριστερό μέλος της εξίσωσης είναι 1. Οπότε δεν χάνουμε ρίζες διαιρώντας την εξίσωσή μας με $\cos x$. Ας δούμε τώρα μαζί τα γραφήματα των συναρτήσεων $y=-\tan x$, και της ευθείας $y=x$.
ef = plot(-tan(x), (x,-3*pi+pi/2, 3*pi-pi/2), ymin=-8, ymax=8, detect_poles = True ,thickness=2);
line = plot(x , (x,-3*pi+pi/2, 3*pi-pi/2), ymin=-8, ymax=8, color = 'green' ,thickness=2)
(ef+line).show(figsize=4)
Από το σχήμα παρατηρούμε ότι σε κάθε υποδιάστημα του πεδίου ορισμού της $y=-\tan x$, δηλαδή σε κάθε ανοιχτό διάστημα της μορφής $(-\frac{\pi}{2}+k\,\pi\,,\,\frac{\pi}{2}+k\,\pi)$, το γράφημα της $y=x$ τέμνει το γράφημα της $y=-\tan x$, σε ένα και μόνο ένα σημείο. Το θέμα είναι τώρα να εντοπίσουμε τα σημεία αυτά, αριθμητικά, αφού συμβολικά αυτό δεν είναι εφικτό. Η εντολή για την αριθμητική επίλυση μιας εξίσωσης στο Sage είναι η find_root. Ας δούμε ποιά συστήματα χρησιμοποιεί το Sage για την αριθμητική επίλυση εξισώσεων.
from sage.misc.citation import get_systems
print get_systems('find_root(x, -1 , 1)')
Το Sage για την αριθμητική εύρεση των ριζών μιας συνάρτησης χρησιμοποιεί το ΝumPy και το SciPy. Πράγματι, εκτελώντας την εντολή find_root?? (με διπλό ερωτηματικό!) μας εμφανίζεται ο κώδικας της εντολής find_root, ο οποίος στο τέλος καταλήγει να καλεί έναν αλγόριθμο του SciPy
import scipy.optimize
return scipy.optimize.brentq(f, a, b,
full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
Οπότε το συμπέρασμα είναι ότι ο αλγόριθμος που χρησιμοποιείται για την αριθμητική εύρεση των ριζών μιας συνάρτησης στο Sage, είναι η κλασική μέθοδο του Brent. Ο αλγόριθμος αυτός συνδυάζει τις μεθόδους της διχοτόμησης, της παρεμβολής και της τετραγωνικής παρεμβολής για τον υπολογισμό της αντίστροφης συνάρτησης της οποίας αναζητούμε τις ρίζες της. Ο αλγόριθμος του Brent είναι μια πολύ δημοφιλής μέθοδος για την εύρεση των ριζών μιας συνάρτησης επειδή είναι ευσταθής, αξιόπιστος και ταχύτατος. Ας δούμε όμως περισσότερες λεπτομέρειες για τον αλγόριθμο του Brent, καλώντας βοήθεια από το SciPy για τον αλγόριθμο
from scipy.optimize import brentq
help(brentq)
Help on function brentq in module scipy.optimize.zeros:
brentq(f, a, b, args=(), xtol=1e-12, rtol=4.4408920985006262e-16, maxiter=100, full_output=False, disp=True)
Find a root of a function in given interval.
.
.
Notes
-----
`f` must be continuous. f(a) and f(b) must have opposite signs.
Διαβάζοντας προσεχτικά παρατηρούμε ότι για να χρησιμοποιήσουμε την μέθοδο του Brent με αξιοπιστία θα πρέπει:
Οπότε τίθεται το ερώτημα, ποιά $f(x)$ θα πάρουμε για να εξασφαλίσουμε την σύγκλιση της μεθόδου αριθμητικά; Μετά από λίγη σκέψη, θα πάρουμε την αρχική μας συνάρτηση $f(x)=\sin x + x\, \cos x$, η οποία είναι συνεχής σε όλο το πεδίο οριμού της δηλαδή το $\mathbb{R}$, και ως διαστήματα θα πάρουμε διαστήματα της μορφής $[-\frac{\pi}{2}+k\,\pi\,,\,\frac{\pi}{2}+k\,\pi]$, στα οποία γνωρίζουμε ότι σε κάθε ένα από αυτά υπάρχει μία και μόνο μία ρίζα της εξίσωσης $f(x)=0$.
Είναι όμως ετερόσημες οι τιμές της συνάρτησης $\,\,f(x)=\sin x+x\,\cos x$, στα άκρα των παραπάνω διαστημάτων $[a,b]$;
f = sin(x)+x*cos(x)
var('k')
assume(k,'odd');
print str('αν k άρτιος:');
print str('f(a)='), f(x=-pi/2 + k*pi).trig_expand() , str('f(b)='), f(x=pi/2 + k*pi).trig_expand()
forget(); assume(k,'even');
print str('αν k περιττός:');
print str('f(a)='), f(x=-pi/2 + k*pi).trig_expand() , str('f(b)='), f(x=pi/2 + k*pi).trig_expand()
forget()
Οπότε, εξασφαλίζονται όλες οι προϋποθέσεις για την επιτυχή εφαρμογή του αλγόριθμου Brent για να βρούμε τις ρίζες της $f(x)$. Με την παρακάτω εντολή βρίσκουμε τις πέντε ρίζες στα αντίστοιχα πέντε διαστήματα της μορφής $[-\frac{\pi}{2}+k\,\pi\,,\,\frac{\pi}{2}+k\,\pi]$ με $k=-2,-1,0,1,2$, και τα τοποθετούμε σε μια λίστα με όνομα w
w=[]
for k in [-2..2]:
root = find_root(f, -pi/2 + k*pi , pi/2 + k*pi)
w.append(root)
print w
Προφανώς, με ανάλογο τρόπο μπορούμε να εντοπίσουμε και τις ρίζες που βρίσκονται στα υπόλοιπα διαστήματα. Για παράδειγμα, ας πάμε στο διάστημα $[-\frac{\pi}{2}+k\,\pi\,,\,\frac{\pi}{2}+k\,\pi]$, με $k=10000$.
j = 10000
root = find_root(f, -pi/2 + j*pi , pi/2 + j*pi,full_output=True)
print root
Tέλος, ας απεικονίσουμε γραφικά τα αποτελέσματά μας στην αρχική γεωμετρική εικόνα των ριζών ως σημεία τομής του γραφήματος της $y=-\tan x$ με την ευθεία $y=x$.
pp=list([ ( w[i] , w[i] ) for i in [0..(len(w)-1)]])
pp_plot=scatter_plot(pp, facecolor='red', marker='o', markersize=30)
(ef+line+pp_plot).show(figsize=4)
Αν εφαρμόζαμε την μέθοδο αριθμητικής εύρεσης των ριζών στην συνάρτηση $g(x)=-\tan x-x$, ο αλγόριθμος θα αποτύγχανε. Ο λόγος είναι ότι δεν υπάρχουν οι τιμές της $g(x)$ στα άκρα των διαστημάτων $[-\frac{\pi}{2}+k\,\pi\,,\,\frac{\pi}{2}+k\,\pi]$, αφού η $g(x)$ δεν ορίζεται στα σημεία αυτά. Επιπλέον, δεν μπορούμε να αποφανθούμε ούτε για το πρόσημο της $g(x)$ στα άκρα των διαστημάτων αφού ούτε τα όρια υπάρχουν εκεί. Για παράδειγμα, στο $x=\frac{\pi}{2}$ τα πλευρικά όρια της $g(x)$ είναι
g = -tan(x) - x;
print g.limit(x=pi/2,dir = '+');
print g.limit(x=pi/2,dir = '-');
brentq(g, pi/2, 3*pi/2 , full_output=true)
6.2.2 Συστήματα υπερβατικών εξισώσεων¶
Προς το παρόν στο Sage δεν έχει υλοποιηθεί ακόμη (version 7.5.1) κάποια εντολή για την αριθμητική επίλυση συστημάτων υπερβατικών εξισώσεων. Όμως, για τον σκοπό αυτό μπορούμε να χρησιμοποιήσουμε αλγόριθμους που είναι διαθέσιμοι στο SciPy σε συνδυασμό με το NumPy. Ορισμένοι αλγόριθμοι χρησιμοποιούν τον Ιακωβιανό πίνακα των συναρτήσεων που ορίζουν το σύστημα που θέλουμε να επιλύσουμε οπότε πριν καλέσουμε τις εντολές από τα SciPy, NumPy μπορούμε να υπολογίσουμε τον πίνακα αυτό συμβολικά από το Sage και να απεικονίσουμε τουλάχιστον γραφικά την λύση των εξισώσεων.
Για παράδειγμα, ας υποθέσουμε ότι θέλουμε να επιλύσουμε το εξής σύστημα εξισώσεων
$$ f(x,y) := e^x + y-4=0\,, \qquad g(x,y) := e^y - x =0 \,, \qquad x,y\in \mathbb{R}\,. $$
Ας ορίσουμε τις συναρτήσεις αυτές στο Sage
var('x y')
f = exp(x)+y-4;
g = exp(y)-x;
print f; print g
Για να απεικονίσουμε γραφικά το(α) σημείο(α) (αν υπάρχουν) που τέμνονται οι καμπύλες $f(x,y)=0$ και $g(x,y)=0$, χρησιμοποιούμε την εντολή contour_plot και από όλες τις ισοσταθμικές καμπύλες $f(x,y)=c_1$ και $g(x,y)=c_2$ επιλέγουμε εκείνες για τις οποίες ισχύει $c_1=c_2=0$.
cf = contour_plot(f, (x,-3,3), (y,-3,3), contours=[0], fill=False, labels=True,cmap='jet', label_inline=True,
label_fmt={0:" f=0 "}, label_colors='black',linewidths=2)
cg = contour_plot(g, (x,-3,3), (y,-3,3), contours=[0], fill=False, labels=True,cmap='autumn', label_inline=True,
label_fmt={0:" g=0 "}, label_colors='black',linewidths=2)
(cf+cg).show(figsize=5)
Παρατηρούμε ότι στο τετράγωνο $S=[-3,3]\times[-3,3]$ υπάρχει ένα σημείο τομής των καμπυλών $f(x,y)=0$, $g(x,y)=0$ και μάλιστα το σημείο αυτό είναι μοναδικό στο $S$. Οπότε, αν ο αλγόριθμος που θα χρησιμοποιήσουμε συγκλίνει και πάρουμε αρχική προσέγγιση κοντά στην μοναδική ρίζα, ο αλγόριθμος μάλλον θα συγκλίνει στην μοναδική ρίζα στο $S$.
Ορισμένοι αλγόριθμοι του SciPy χρησιμοποιούν τον Ιακωβιανό πίνακα $J(f,g;x,y)$, οπότε το Sage μπορεί να μας βοηθήσει να υπολογίσουμε συμβολικά τον πίνακα $J(f,g;x,y)$ και να τον περάσουμε ακριβώς στους αλγόριθμους αυτούς, αντί ο αντίστοιχος αλγόριθμος του SciPy να τον υπολογίσει προσεγγιστικά. Επίσης, επειδή το SciPy χρησιμοποιεί την Python, πρέπει να μετατρέψουμε τις συναρτήσεις $(f,g)$ και τον πίνακα $J$ σε ανάλογο συμβολισμό. Αυτά υλοποιούνται στις παρακάτω εντολές οι οποίες καταλήγουν στην εύρεση της αριθμητικής τιμής της ρίζας που ψάχνουμε.
F = (f,g)
print jacobian(F,(x,y))
import numpy as np
from scipy import optimize
def fun(x):
return [ exp(x[0]) + x[1] - 4.0, exp(x[1])-x[0] ]
def jac(x):
return np.array([[ exp(x[0]) , 1.0],
[ -1.0 , exp(x[1])]])
sol = optimize.root(fun, [1, 1] , jac = jac , method='hybr'); print sol.x
Η αρχική εκτίμηση της ρίζας που δώσαμε ήταν το σημείο $(1,1)$, και ο αλγόριθμος που χρησιμοποιήσαμε είναι η hybr (μια τροποποιημένη μεθόδος του αλγόριθμου του Powell όπως μας πληροφορεί στην σχετική βοήθεια το SciPy). Την αριθμητική τιμή της λύσης που μας έδωσε το SciPy μπορούμε να την καταχωρίσουμε στο Sage για περαιτέρω μελέτη. Για παράδειγμα
x0 = sol.x[0]; y0 = sol.x[1]; print str('x0='), x0 ; print str('y0='), y0 ;
print f(x=x0,y=y0);
print g(x=x0,y=y0),
από όπου συμπεραίνουμε ότι η λύση που βρήκαμε ικανοποιεί και τις δυο εξισώσεις $f(x,y)=0$, $g(x,y)=0$, σχεδόν στην ακρίβεια δεκαδικών ψηφίων κινητής υποδιαστολής που εκτελεί ο υπολογιστής μας. Για διαφορετικούς αλγόριθμους αριθμητικής επίλυσης συστημάτων εξισώσεων με το SciPy παραπέμπουμε τον αναγνώστη στην σχετική βοήθεια.
help(optimize.root)
6.3 Συστήματα πολυωνυμικών εξισώσεων
6.3.1 Βάσεις Gröbner¶
Αν το σύστημα των εξισώσεων που θέλουμε να επιλύσουμε απαρτίζεται από πολυώνυμα τότε το Sage μας δίνει την δυνατότητα να βρούμε μια βάση Gröbner. Η κατασκευή μιας βάσης Gröbner για ένα πολυωνυμικό ιδεώδες είναι το μη-γραμμικό ανάλογο της διαδικασίας να φέρουμε ένα σύστημα γραμμικών εξισώσεων σε ανηγμένη κλιμακωτή μορφή. Η ύπαρξη μιας βάσης Gröbner μπορεί να συναχθεί από το λεγόμενο Θεώρημα βάσης του Hilbert. Όμως η εύρεση με ρητό τρόπο μιας βάσης Gröbner μπορεί να είναι ένα αρκετά δύσκολο κι επίπονο εγχείρημα. Ο λεγόμενος αλγόριθμος του Buchberger δίνει την γενική λύση στην κατασκευή μιας βάσης Gröbner από ένα ιδεώδες πολυωνύμων κι έχει υλοποιηθεί στο Sage μέσω της Singular.
Μια συστηματική διαπραγμάτευση των βάσεων Gröbner μας βγάζει πολύ έξω από τους στόχους των σημειώσεων αυτών. Αντί αυτής της διαπραγμάτευσης θα παρουσιάσουμε ένα παράδειγμα από την γεωμετρία που μας δίνει μια εποπτική εικόνα επίλυσης πολυωνυμικών εξισώσεων μέσω της κατασκευής μιας βάσης Gröbner.
Το θεώρημα του Descartes (1643) λέει ότι οι ακτίνες οποιωνδήποτε τεσσάρων κύκλων οι οποίοι εφάπτονται αμοιβαία και οι τέσσερις, ικανοποιεί μια συγκεκριμένη δευτεροβάθμια εξίσωση. Δοσμένων τριών αμοιβαία εφαπτόμενων κύκλων, λύνοντας την δευτεροβάθμια εξίσωση μπορούμε να κατασκευάσουμε τον τέταρτο κύκλο.
Είναι γνωστό από την Αναλυτική Γεωμετρία ότι δυο κύκλοι $(O_1(x_1,y_1),r_1)$ και $(Ο_2(x_2,y_2),r_2)$ με κέντρα $O_1$, $O_2$ και ακτίνες $r_1,r_2$, αντίστοιχα, που εφάπτονται εξωτερικά ικανοποιούν την εξίσωση
$$ r_1+r_2 = {\overline{O_1O_2}} \quad \Rightarrow \quad r_1+r_2 = \sqrt{(x_1-x_2)^2 + (y_1-y_2^2)}\,.$$
Επιλέγουμε τρια τυχαία σημεία στο επίπεδο εκ των οποίων τα δύο είναι πάνω στο x-άξονα. Επειδή θέλουμε οι κύκλοι που θα έχουν κέντρα τα σημεία αυτά να εφάπτονται εξωτερικά επιλέγουμε τα σημεία αυτά ανάλογα.
(c1x , c1y) = (RR.random_element(0.5, 1).n(digits=2) , 0)
(c2x , c2y) = (RR.random_element(-1, -0.5).n(digits=2) , 0)
(c3x , c3y) = (RR.random_element(-1, 1).n(digits=2) , RR.random_element(1, 1.5).n(digits=2))
p1 = point((c1x, c1y), size=25, color='blue')
p2 = point((c2x, c2y), size=25, color='blue')
p3 = point((c3x, c3y), size=25, color='blue')
(p1+p2+p3).show(figsize=3,axes=False)
Στην συνεχεία καταστρώνουμε τις εξισώσεις που μας δίνουν τις ακτίνες $r_1,r_2,r_3$ των τριών κύκλων, υποθέτοντας ότι αυτοί εφάπτονται αμοιβαία.
var('r1 , r2, r3')
eq1 = r1 + r2 == sqrt((c1x-c2x)^2 + (c1y-c2y)^2);
eq2 = r1 + r3 == sqrt((c1x-c3x)^2 + (c1y-c3y)^2);
eq3 = r2 + r3 == sqrt((c2x-c3x)^2 + (c2y-c3y)^2);
sol=solve([eq1,eq2,eq3],[r1,r2,r3],solution_dict=true)
Το παραπάνω σύστημα είναι ένα γραμμικό σύστημα για τις ακτίνες το οποίο λύνεται εύκολα και καταχωρούμε τις λύσεις στις μεταβλητές των ακτίνων
rr1 = sol[0][r1];
rr2 = sol[0][r2];
rr3 = sol[0][r3];
Το αποτέλεσμα είναι οι εξής τρεις αμοιβαία εφαπτόμενοι κύκλοι
var('x,y')
c1 = implicit_plot((x - c1x)^2 + (y - c1y)^2 - rr1^2, (x, -3, 3), (y, -3, 3), color='green')
c2 = implicit_plot((x - c2x)^2 + (y - c2y)^2 - rr2^2, (x, -3, 3), (y, -3, 3), color='green')
c3 = implicit_plot((x - c3x)^2 + (y - c3y)^2 - rr3^2, (x, -3, 3), (y, -3, 3), color='green')
(p1 + p2 + p3 + c1 + c2 + c3 ).show(figsize=6,frame=False)
Το πρόβλημα του Descartes, το οποίο αποτελεί ένα περιορισμένο πρόβλημα του Απολλώνιου από την Πέργαμο, είναι να βρούμε έναν τέταρτο κύκλο ο οποίος είναι εφαπτόμενος και στους τρεις δοσμένους αμοιβαία εφαπτόμενους κύκλους.
Ας ονομάσουμε με $Ο_4(x,y)$ και $r$, το κέντρο και την ακτίνα, αντίστοιχα, του τέταρτου κύκλου. Τότε είναι προφανές ότι για να εφάπτεται ο τέταρτος κύκλος αμοιβαία και στους τρεις άλλους κύκλους θα πρέπει να ισχύει το σύστημα των εξισώσεων
$$(x - x_1)^2 + (y - y_1)^2 - (r + r_1)^2 =0 \,,$$
$$(x - x_2)^2 + (y - y_2)^2 - (r + r_2)^2 =0 \,,$$
$$(x - x_3)^2 + (y - y_3)^2 - (r + r_3)^2 =0 \,.$$
Το παραπάνω πρόκειται για ένα σύστημα τριών πολυωνυμικών εξισώσεων δευτέρου βαθμού για τις τρεις άγνωστες ποσότητες, το κέντρο $(x,y)$ και την ακτίνα $r$ του τέταρτου κύκλου $O_4$. Θα επιλύσουμε το σύστημα αυτό στο Sage κατασκευάζοντας μια βάση Gröbner. Πρώτα απ' όλα σχηματίζουμε συμβολικά το σύστημα των εξισώσεων.
var('r')
e1 = (x - c1x)^2 + (y - c1y)^2 - (r + rr1)^2
e2 = (x - c2x)^2 + (y - c2y)^2 - (r + rr2)^2
e3 = (x - c3x)^2 + (y - c3y)^2 - (r + rr3)^2
eqs = [e1, e2, e3]
Στην συνέχεια δημιουργούμε έναν δακτύλιο πολυωνύμων τριών μεταβλητών $(x,y,r)$ με συντελεστές στους ρητούς, με λεξικογραφική διάταξη μονονύμων, και μεταγράφουμε στον δακτύλιο αυτό, τα πολυώνυμα που ορίζουν τις εξισώσεις μας.
PR.<x,y,r> = PolynomialRing(QQ, order='lex')
print PR
q1 = PR(e1); q2 = PR(e2); q3 = PR(e3)
print q1;
Κατόπιν, τα πολυώνυμα $q_1,q_2,q_3$ τα θέτουμε σε μια λίστα και ορίζουμε το ιδεώδες των πολυωνύμων για το οποίο θα κατασκευάσουμε μια βάση Gröbner
qs = [q1, q2, q3]
I = Ideal(qs)
print I
Με την εντολή groebner_basis() το Sage μας επιστρέφει την βάση Gröbner του $Ι$
B = I.groebner_basis()
print B
Παρατηρούμε ότι το τρίτο πολυώνυμο στην λίστα (βάση Gröbner) που μας έδωσε το Sage είναι δευτέρου βαθμού και αφορά μόνο την μεταβλητή $r$, την ακτίνα του τέταρτου κύκλου. Τα πρώτα δυο πολυώνυμα είναι γραμμικά στις μεταβλητές $x,y,r$. Οπότε λύνοντας την δευτεροβάθμια εξίσωση για το $r$ παίρνουμε δυο λύσεις (ένα εσωτερικό κι έναν εξωτερικό εφαπτόμενο κύκλο) και αντικαταθιστώντας στις δυο πρώτες παίρνουμε τα αντίστοιχα κέντρα των κύκλων. Πραγματικά
soln = I.variety(RR)
(x1s , y1s , r1s) = (soln[0][x] , soln[0][y] , soln[0][r])
(x2s , y2s , r2s) = (soln[1][x] , soln[1][y] , soln[1][r])
var('x , y')
kyk1=implicit_plot((x-x1s)^2+(y-y1s)^2-r1s^2, (x, -3, 3), (y, -3, 3), color='red')
kyk2=implicit_plot((x-x2s)^2+(y-y2s)^2-r2s^2, (x, -3, 3), (y, -3, 3), color='red')
(p1 + p2 + p3 + c1 + c2 + c3 + kyk1 + kyk2).show(figsize=7,frame=False)
6.3.2 Απαλείφουσα¶
Μια εναλλακτική μέθοδος για την επίλυση πολυωνυμικών εξισώσεων είναι να χρησιμοποιήσουμε την απαλείφουσα, η οποία είναι διαθέσιμη στο Sage μέσω της Singular. Η απαλείφουσα είναι μια ορίζουσα και αποτελεί το μη-γραμμικό ανάλογο της γνωστής μας ορίζουσας για γραμμικά συστήματα εξισώσεων. Πιο συγκεκριμένα.
Η απαλείφουσα δυο πολυωνύμων
$$ f = a_m\, x^m + a_{m-1}\,x^{m-1} + \cdots a_0\,,$$
$$ g = b_n\, x^n + b_{n-1}\,x^{n-1} + \cdots b_0\,,$$
με $a_m,b_n\neq 0 $, ορίζεται ως η $(m+n)\times (m+n)$ ορίζουσα
$$ {\mathrm{Res}}(f,g) = \det
\left(\begin{array}{cccccccc}
a_m & a_{m-1} & \cdots & a_0 & 0 & 0 & \ldots & 0 \\
0 & a_{m} & a_{m-1} & \cdots & a_0 & 0 & \ldots & 0 \\
0 & 0 & \ddots & \ddots & \cdots & \vdots & \vdots & 0 \\
0 & 0 & \cdots & 0 & a_m & a_{m-1} & \ldots & a_0 \\
b_n & b_{n-1} & \cdots & b_0 & 0 & 0 & \ldots & 0 \\
0 & b_{n} & b_{n-1} & \cdots & b_0 & 0 & \ldots & 0 \\
0 & 0 & \ddots & \ddots & \cdots & \vdots & \vdots & 0 \\
0 & 0 & \cdots & 0 & b_n & b_{n-1} & \cdots & b_0 \\
\end{array} \right)$$
Οι πρώτες $n$-γραμμές αφορούν τους συντελεστές $a_i$, και οι τελευταίες $m$-γραμμές αφορούν τα $b_i$. Καλείται ορισμένες φορές και ως απαλείφουσα του Sylvester.
Για παράδειγμα, στο Sage η απαλείφουσα των πολυωνύμων
$$ f(x) = a\,x^2 + b\,x+c\,,\qquad g(x)=d\,x+e\,, $$
υπολογίζεται από την ορίζουσα του παραπάνω πίνακα Sylvester, ο οποίος παράγεται χρησιμοποιώντας την Singular για να δηλώσουμε ένα δακτύλιο πολυωνύμων
και με την εντολή .sylvester_matrix(). Συγκεκριμένα,
PR.<x,a,b,c,d,e> = PolynomialRing(QQ)
print PR
f = a*x^2 + b*x + c
g = d*x + e
print type(f), type(g)
Παρατηρούμε ότι και τα δυο πολυώνυμα $f$ και $g$ έχουν τύπο MPolynomial_libsingular.
Res_x = f.sylvester_matrix(g, x)
print Res_x
print Res_x.det()
Τα βασικά θεωρήματα σχετικά για την απαλείφουσα είναι:
Για παράδειγμα, όπως διαπιστώσαμε πιο πάνω η απαλείφουσα των πολυωνύμων $ f(x) = a\,x^2 + b\,x+c$, και $g(x)=d\,x+e\,, $ είναι $a\,e^2 - b\,d\,e + c\,d^2 $. Το τελευταίο είναι ένα πολυώνυμο πέντε μεταβλητών στον δακτύλιο $k[a,b,c,d,e]$, με συντελεστές σε οποιοδήποτε σώμα $K$. Αυτό που λέει το θεώρημα είναι ότι το πολυώνυμο $a\,e^2 - b\,d\,e + c\,d^2 $ είναι ανάγωγο στο $Κ$ (επομένως και στο $k[a,b,c,d,e]$).
Ας δούμε πως μπορούμε να επιλύσουμε ένα σύστημα πολυωνυμικών εξισώσεων με το Sage στην πράξη.
Θεωρούμε τα εξής πολυώνυμα σε δυο μεταβλητές $x,y$:
$$p(x,y) = x^2 + y^2 - 5 \,, \qquad q(x,y) = x^2 - y^2 - 3 \,.$$
Θέλουμε να βρούμε τις λύσεις του συστήματος πολυωνυμικών εξισώσεων $p(x,y)=0$, $q(x,y)=0$.
Θα απαλείψουμε την μεταβλητή $x$ από το πολυώνυμο p, ως προς το πολυώνυμο q, ή ισοδύναμα θα απαλείψουμε την μεταβλητή $x$ από το πολυώνυμο q, ως προς το πολυώνυμο p.
PR.<x,y,w> = PolynomialRing(QQ)
print PR
p = x^2 + y^2 - 5
q = x^2 - y^2 - 3
print type(p), type(q)
print p.resultant(q,x)
print q.resultant(p,x).factor()
Ο μηδενισμός της απαλείφουσας σημαίνει ότι θα πρέπει $y=1$ ή $y=-1$. Θέτουμε τις τιμές αυτές στο $p$ (ή ισοδύναμα στο $q$) και λύνουμε το πολυώνυμο που προκύπτει, πλέον ως προς μία μεταβλητή την $x$.
R = QQ[x]
p1 = R(p(y=1)); print p1;
print p1.roots(ring=CC)
R = QQ[x]
p_red = R(p(y=-1)); print p_red;
print p_red.roots(ring=CC)
Από τα παραπάνω συμπεραίνουμε ότι οι λύσεις του συστήματος πολυωνυμικών εξισώσεων $p(x,y)=0$, $q(x,y)=0$, είναι οι εξής: $$ (-2,1)\,,\quad (2,1)\,\quad (-2,-1)\,,\quad (2,-1)\,.$$
Χρησιμοποιώντας με παρόμοιο τρόπο την απαλείφουσα σε συστήματα πολυωνυμικών εξισώσεων με περισσότερες μεταβλητές μπορούμε να απαλείφουμε διαδοχικά μεταβλητές, με στόχο να πάρουμε μια εξίσωση μιας μεταβλητής και αφού λύσουμε την τελευταία να επιλύσουμε ως προς τις υπόλοιπες μεταβλητές με την ίδια διαδικασία.
6.4 Επίλυση ανισώσεων - Βελτιστοποίηση¶
6.4.1 Συμβολική και αριθμητική επίλυση ανισώσεων - Γραμμικός προγραμματισμός¶
Το Sage μπορεί να λύσει συμβολικά μέσω του Maxima και συστήματα ανισώσεων. Για παράδειγμα
var('x y')
print solve_ineq( x^2 - x > 5)
print solve_ineq([x-y<0,x+y-3<0],[y,x])
print solve_ineq([x-y<0,x+y-3<0],[x,y])
Όταν δεν εφικτό να έχουμε αποτέλεσμα σε συμβολική μορφή, το Sage θα προσπαθήσει να επιλύσει την ανίσωση αριθμητικά.
print solve_ineq(-20 * x - 30 <= 4 * x^3 - 7 * x, x)
Οι μέθοδοι βελτιστοποίησης μιας συνάρτησης έχουν πολύ ενδιαφέρουσες εφαρμογές σε προβλήματα των μαθηματικών και της μηχανικής όπου χρησιμοποιούνται στην προσαρμογή της συνάρτησης σε δοσμένα σύνολα και στην εύρεση της βέλτιστης (μέγιστης/ελάχιστης) τιμής της ως προς τις σχεδιαστικές παραμέτρους. Γενικά η βελτιστοποίηση μη-γραμμικών προβλημάτων είναι αρκετά δύσκολο πρόβλημα. Ας δούμε πως μπορούμε να επιλύσουμε με την βοήθεια του Sage ένα απλό βαθμωτό γραμμικό πρόβλημα βελτιστοποίησης όπου οι παράμετροι παίρνουν τιμές στους πραγματικούς αριθμούς.
Πιο συγκεκριμένα θέλουμε να βρούμε την λύση του ακόλουθου προβλήματος:
Να βρεθεί το
$$ \max \,(x+5\,y)$$
με τον περιορισμό ότι
$$
\begin{cases} x + 0.2 \, y & \leq 4 \,,\\
1.5 \, x+3\,y & \leq 4 \,.
\end{cases} $$
όπου $x$, $y$ είναι πραγματικοί θετικοί αριθμοί. Στο Sage το πρόβλημα αυτό επιλύεται ως εξής:
p = MixedIntegerLinearProgram(maximization=True)
x = p.new_variable(nonnegative=True)
p.set_objective(x[1] + 5*x[2])
p.add_constraint(x[1] + 0.2*x[2], max=4)
p.add_constraint(1.5*x[1] + 3*x[2], max=4)
print p.solve()
από όπου συμπεραίνουμε ότι η μέγιστη τιμή που ψάχνουμε είναι 6.66666666667. Για να δούμε με περισσότερες λεπτομέρειες πως το Sage μπορεί να μας βοηθήσει να επιλύσουμε ένα πρόβλημα γραμμικού προγραμματισμού. Πρώτα απ' όλα μπορεί να μας δώσει μια γραφική αναπαράσταση της λύσης που ψάχνουμε.
A = ([1, 0.2], [1.5, 3])
b = (4, 4)
c = (1, 5)
P = InteractiveLPProblem(A, b, c, ["x1", "x2"], problem_type="max", \
constraint_type="<=", variable_type=">=")
pp = P.plot(); pp.show(figsize=6)
Το Sage χώρισε το επίπεδο σε περιοχές ανάλογα με τις ανισότητες των περιορισμών του προβλήματος, οι οποίες συναληθεύουν στην τριγωνική περιοχή που συμβολίζει με $F$. Επιπλέον σχεδίασε ευθείες παράλληλες (διακεκομμένες) με κλίση $-1/5$, όπου πάνω στις ευθείες αυτές η τιμή της συνάρτησης $x+5\,y$ που θέλουμε να μεγιστοποιήσουμε παραμένουν σταθερές και ίσες με $c$, $x+5\,y=c$.
Όσο απομακρυνόμαστε από την αρχή των αξόνων προς τα πάνω οι τιμές της $c$ μεγαλώνουν, οπότε το σημείο τομής των ευθειών $x+5\,y=c$ (ευθεία με έντονο μαύρο) με τις ευθείες $1.5\,x+3\,y=4$ (ευθεία με πράσινο), και $x=0$ (ευθεία με τυρκουάζ), είναι το ζητούμενο σημείο στο οποίο η $x+5\,y$ παίρνει μέγιστη στην περιοχή $F$. Πραγματικά, από το σχήμα που μας έδωσε το Sage μπορούμε αμέσως να βρούμε ακριβώς το σημείο τομής και την μέγιστη τιμή που ψάχνουμε
var('x y c')
sol = solve([x+5*y==c,1.5*x+3*y==4,x==0],x,y,c,solution_dict=True); print sol
print sol[0][c].n()
Η μέγιστη τιμή της συνάρτησης $x+5\,y$ είναι $c=20/3$, στο σημείο $(x,y)=(0,4/3)$. Το Sage μπορεί να επιλύσει το πρόβλημα αριθμητικά χρησιμοποιώντας την μέθοδο simplex το οποίο είναι πολύ χρήσιμο όταν έχουμε να αντιμετωπίσουμε γραμμικά προβλήματα πολλών μεταβλητών. Με την παρακάτω εντολή μπορούμε να δούμε το αποτέλεσμα κάθε βήματος της μεθόδου simplex στην μορφή κειμένου και εξισώσεων $\mathrm{\LaTeX}$
P = P.standard_form()
P.run_simplex_method()
6.4.2 Αδέσμευτα ακρότατα¶
Ολοκληρώνουμε το κεφάλαιο αυτό παρουσιάζοντας την επίλυση ενός προβλήματος εύρεσης των ακροτάτων μιας συνάρτησης δυο μεταβλητών, χρησιμοποιώντας το Sage, τόσο συμβολικά όσο και αριθμητικά.
Θεωρούμε την συνάρτηση $f:\mathbb{R}^2 \rightarrow \mathbb{R}$ με τύπο
$$ f(x,y) = (3+x-y^2)^2 + (x - 1)^2 + (y - 1)^2 \,.$$
Στο Sage ορίζουμε την συνάρτηση $f(x,y)$ και κάνουμε μια γραφική παράστασή της
var('x,y')
f(x,y) = (3+x-y^2)^2 + (x - 1)^2 + (y - 1)^2
pf = plot3d(f(x,y), (x, -2, 3), (y, -2, 3), adaptive=True, rgbcolor=(1,0.7,0), opacity=0.5)
Για να βρούμε τα υποψήφια τοπικά ακρότατα της $f$ στο $\mathbb{R}^2$ θα πρέπει να επιλύσουμε το σύστημα των εξισώσεων που προκύπτει από τον μηδενισμό της κλίσης $\nabla f = (f_x,f_y)$. Έχουμε
fx=f.diff(x); fy=f.diff(y);
eq1 = fx==0;
eq2 = fy==0;
print fx; print fy
akro = solve([eq1, eq2],x,y,solution_dict=true) ; print akro
Υπάρχουν τρία υποψήφια τοπικά ακρότατα. Για το πρώτο σημείο στην λίστα akro έχουμε ότι η Hessian στο σημείο είναι αρνητική
D = det(f.diff().diff())
print D(x=akro[0][x],y=akro[0][y])
Άρα από γνωστό μας κριτήριο το σημείο αυτό δεν είναι ούτε μέγιστο ούτε ελάχιστο, είναι σαγματικό. Για το δεύτερο και τρίτο σημείο έχουμε ότι οι αντίστοιχες τιμές της Hessian και η δεύτερη παράγωγος $f_{xx}$ στα σημεία αυτά είναι
print D(x=akro[1][x],y=akro[1][y])
print f.diff(x,2)(x=akro[1][x],y=akro[1][y])
print D(x=akro[2][x],y=akro[2][y])
print f.diff(x,2)(x=akro[2][x],y=akro[2][y])
Οπότε από το κριτήριο χαρακτηρισμού τοπικών ακροτάτων τα σημεία αυτά είναι τοπικά ελάχιστα. Απεικονίζουμε και τα τρία σημεία πάνω στην επιφάνεια για να έχουμε μια καλύτερη εποπτεία που ακριβώς βρίσκονται. Με κόκκινο χρώμα απεικονίζουμε το σαγματικό σημείο ενώ με μαύρο τα τοπικά ελάχιστα.
Εναλλακτικά, μπορούμε να λύσουμε το σύστημα των εξισώσεων κατασκευάζοντας μια βάση Gröbner, ως εξής
PR.<x,y> = PolynomialRing(QQ, order='lex')
print PR
q1 = PR(eq1); q2 = PR(eq2);
print q1; print q2;
qs = [q1, q2]
I = Ideal(qs)
print I
B = I.groebner_basis()
print B
R = RealField(30*log(10,2))
akroGroebner = I.variety(R); print akroGroebner
όπου χρησιμοποιήσαμε ένα σώμα πραγματικών με μεγάλη ακρίβεια.
Ας υποθέσουμε τώρα ότι δεν ήταν εφικτή η συμβολική επίλυση του συστήματος των εξισώσεων, όμως θέλαμε να έχουμε μια αριθμητική εκτίμηση των τοπικών ελαχίστων. Για τέτοιες περιπτώσεις το Sage μας παρέχει μέσω της εντολής minimize έναν αλγόριθμο που εκτελεί αριθμητική εύρεση των ελαχίστων της $f(x,y)$, όπου όμως θα πρέπει να το τροφοδοτήσουμε με μια αρχική εκτίμηση. Για το συγκεκριμένο παράδειγμα έχουμε
p = minimize(f, x0=[-1, -1])
print 'minimum at', p, 'with value', f(p[0], p[1])
p = minimize(f, x0=[1, 1])
print 'minimum at', p, 'with value', f(p[0], p[1])
Παρατηρούμε ότι οι αριθμητικές προσεγγίσεις των ελαχίστων συμφωνούν με τα ακριβή αποτελέσματα που πήραμε με την μέθοδο των βάσεων Gröbner στα πρώτα 6 έως 8 δεκαδικά ψηφία.