MEM282: Mathematical Modeling

Πληθυσμιακά μοντέλα

Άσκηση (Ανάπτυξη πληθυσμού)
Θεωρήστε ότι ένας πληθυσμός N ζώων ζουν σε ένα δάσος. Κάθε χρόνο το 5% των ζευγαριών αποκτούν ένα παιδί.
(α) Γράψτε ένα πρόγραμμα το οποίο να υπολογίζει τον αριθμό των ζώων κάθε χρόνο για τα επόμενα n χρόνια. Κάνετε γραφική παράσταση N=N(i), όπου i είναι το έτος.
Ακολούθως θεωρήστε ότι κάθε χρόνο πεθαίνουν 1.5% των ζώων.
(β) Γράψτε ένα πρόγραμμα το οποίο να υπολογίζει τον αριθμό των ζώων κάθε χρόνο για τα επόμενα n χρόνια. Κάνετε γραφική παράσταση.

Άσκηση (Απλό μοντέλο ανάπτυξης πληθυσμού)
Μελετήστε το απλό μοντέλο ανάπτυξης ενός πληθυσμού N(t): \[ \frac{dN}{dt} = -a N, \qquad a >0, \] (α) Λύστε αριθμητικά την εξίσωση για αρχική συνθήκη N(t=0)=N₀ της επιλογής σας.
(β) Κάνετε γραφική παράσταση του πληθυσμού ως συνάρτηση του χρόνου N=N(t). (γ) Ποιά η αναλυτική λύση της εξίσωσης; Κάνετε σύγκριση με την αριθμητική λύση.

(α) Δείτε τον κώδικα matlab

(β,γ) Μπορείτε να βελτιώσετε το παραπάνω πρόγραμμα και την παραγόμενη γραφική παράσταση με τους εξής τρόπους:

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

Δείτε τον βελτιωμένο κώδικα matlab και τον αντίστοιχο κώδικα python.

Η γραφική παράσταση από τον βελτιωμένο κώδικα matlab:

simple growth model

Δείτε τη γραφική παράσταση των αποτελεσμάτων με το plot.ly εδώ.

Ερωτήσεις:

Χρησιμοποιήστε μικρό a (π.χ., a=0.1), δείτε τι αλλαγή επιφέρει στη λύση και προσαρμόστε κατάλληλα το γράφημά σας.
Αντικαταστήστε το δεξιό μέλος της εξίσωσης με $aN$ και δείτε τι αποτέλεσμα προκύπτει (κάνετε γραφική παράσταση).

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

Άσκηση (Μοντέλο λογιστικής ανάπτυξης πληθυσμού)
Μελετήστε το μοντέλο λογιστικής ανάπτυξης ενός πληθυσμού N(t): \[ \frac{dN}{dt} = a \left( 1 - \frac{N}{K} \right) N, \] όπου a και K είναι θετικές σταθερές.
(α) Γράψτε την εξίσωση σε αδιάστατη μορφή.
(β) Λύστε αριθμητικά την εξίσωση για αρχική συνθήκη N(t=0)=N₀ της επιλογής σας.
(γ) Σχεδιάστε τον πληθυσμό ως συνάρτηση του χρόνου N(t) και συγκρίνετε την αριθμητική με την ακριβή λύση.

(α) Ορίζουμε τις μεταβλητές τ=? και n=? και έχουμε
dn/dτ = (1-n) n

(β) Δείτε τον κώδικα matlab εδώ και την απαραίτητη function εδώ.

Παρατηρήσεις:

Ξέρουμε ότι αν επιλέξουμε αρχική συνθήκη N(t=0) > K τότε στη χρονική εξέλιξη το N θα μειωθεί μέχρι να φθάσει στην τιμή K. Αντιστοίχως, αν επιλέξουμε N(t=0) < K στη χρονική εξέλιξη το N θα αυξάνεται μέχρι να φθάσει στην τιμή K.
Ο ρυθμός με τον οποίο μεταβάλλεται το N είναι ανάλογος περίπου του 1/a. Με αυτό το κριτήριο επιλέγουμε το t_fin >> 1/a, ώστε το σύστημα να έχει (σχεδόν) φθάσει στο σημείο ισορροπίας για t=t_fin.

(γ) Η ακριβής λύση είναι: N(t) = [N₀ K exp(at)]/[K + N₀ (exp(at)-1)]
και στις αδιάστατες μεταβλητές είναι: n(τ) = [n₀ exp(τ)]/[1 + n₀ (exp(τ)-1)]

Η γραφική παράσταση από τον κώδικα matlab:

logistic growth model

Ερωτήσεις.

Χρησιμοποιήστε αρχικό πληθυσμό μεγαλύτερο του K και δείτε τι αποτέλεσμα προκύπτει.
Προσπαθήστε να καταλάβετε σε τι διαφέρει το γράφημα με τις αρχικές από αυτό με τις αδιάστατες μεταβλητές.

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

Άσκηση (Μοντέλο Lotka-Volterra)
Μελετήστε το μοντέλο Lotka-Volterra για τους πληθυσμούς δύο ειδών x, y: \begin{align} \dot{x} & = a x - c x y \\ \dot{y} & = -b y + d x y. \end{align} (α) Βρείτε το σημείο ισορροπίας του συστήματος (κάνετε μία δική σας επιλογή για τις σταθερές a, b, c, d > 0).
(β) Βρείτε αριθμητική λύση x(t), y(t) (για αρχική συνθήκη της επιλογής σας) και σχεδιάστε τους πληθυσμούς ως συναρτήσεις του χρόνου x=x(t), y=y(t).
(γ) Βρείτε αριθμητικές λύσεις και σχεδιάστε αρκετές καμπύλες στον χώρο των φάσεων (x, y).
(δ) Σχεδιάστε τις καμπύλες στο χώρο των φάσεων χρησιμοποιώντας το ολοκλήρωμα του συστήματος εξισώσεων.

(α) Θέτουμε \begin{align} a x - c x y = 0 \\ -b y + d x y = 0. \end{align} και βρίσκουμε ως λύση το σημείο $(x_0, y_0) = (b/d, a/c)$.

(β) Μπορείτε να βρείτε τον απαραίτητο κώδικα matlab στον ιστότοπο της Mathworks, εδώ.
Επίσης και στο [BF], σελ. 109.

Ο δικός μας κώδικας matlab βρίσκεται εδώ και η απαραίτητη function εδώ.
Ο δικός μας κώδικας python βρισκεται εδώ.

Προσέξτε τα ακόλουθα σημεία:

Επιλέγοντας τις τιμές των παραμέτρων a,b,c,d καθορίσουμε και τη θέση του σημείου ισορροπίας. Αυτό μπορεί να αποτελέσει κριτήριο για την επιλογή των παραμέτρων.
Οι λύσεις του συστήματος είναι περιοδικές. Ώστε, θα πρέπει να επιλέξουμε το χρονικό διάστημα στο οποίο θα λύσουμε το σύστημα ώστε να είναι μεγαλύτερο της περιόδου των λύσεων που μας ενδιαφέρουν.
Επιλέγουμε τις αρχικές συνθήκες σε σχέση με το σημείο ισορροπίας, το οποίο είναι κέντρο. Αναμένουμε να πάρουμε περιοδικές λύσεις που θα ταλαντώνονται γύρω από το σημείο ισορροπίας.

Ο κώδικας δίνει τις ακόλουθες γραφικές παραστάσεις (για δύο διαφορετικές αρχικές συνθήκες):
Lotka-Volterra solution

(γ) Επαναλαμβάνουμε τον αριθμητικό υπολογισμό για αρκετές αρχικές συνθήκες και σχεδιάζουμε καμπύλες στο επίπεδο (x, y). Η τακτική που ακολουθούμε για την επιλογή αρχικών συνθηκών και έλεγχο των αποτελεσμάτων είναι η εξής:

προετοιμασία: αναφερόμαστε στις αναλυτικές πράξεις που έχουμε κάνει για το συγκεκριμένο μοντέλο,
υπολογισμοί: επιλέγουμε παραμέτρους και υπολογίζουμε πού βρίσκεται το σημείο ισορροπίας και τι είδους σημείο είναι (π.χ., κέντρο),
αρχικές συνθήκες: επιλέγουμε στον κώδικα αρχικές συνθήκες γύρω από το σημείο ισορροπίας,
επιβεβαίωση: τρέχουμε τον κώδικα και πρέπει να βγάλει μία τροχιά γύρω από το σημείο ισορροπίας που να επιβεβαιώνει ότι το σημείο είναι, π.χ., κέντρο (δηλαδή, στο Lotka-Volterra, αναμένουμε να πάρουμε, στο διάγραμμα φάσεων, τροχιές κλειστές γύρω από το σημείο ισορροπίας.

Ο κώδικας βρίσκεται εδώ.

Παίρνουμε το ακόλουθο αποτέλεσμα στον χώρο φάσεων:
Lotka-Volterra phase portrait

(δ) Για να σχεδιάσουμε το διάγραμμα φάσεων μπορούμε επίσης να χρησιμοποιήσουμε το ολοκλήρωμα το οποίο έχουμε βρει για το σύστημα Lotka-Volterra: \[ (c\,y - a \ln y) + (d\,x - b \ln x) = C. \] Πρέπει να σχεδιάσουμε ισοσταθμικές καμπύλες της συνάρτησης $H(x,y) = (c\,y - a \ln y) + (d\,x - b \ln x)$. Αυτό μπορούμε να το κάνουμε με το gnuplot, matlab, plot.ly.

Για παράδειγμα, στο gnuplot δίνουμε τις εντολές:
a=10; c = 4; b = 2; d = 1 # (δίνουμε τιμές στις παραμέτρους)
set contour; unset surface # (θα σχεδιάσουμε ισοσταθμικές και όχι επιφάνεια)
set isosamples 100,10 # (αριθμός σημείων κάθε γραμμής που σχεδιάζεται)
set xrange[0:7]; set yrange[0:7] # (όρια αξόνων)
set view 0,0 # (οπτική γωνία για τη γραφική παράσταση)
et cntrparam levels incremental 2,3,17 # (contour levels: 2,5,8,11,17)
splot c*y-a*log(y) + d*x - b*log(x) # (σχεδιάζουμε τις ισοσταθμικές)

Ερωτήσεις.

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

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

Άσκηση (Μοντέλο ανταγωνισμού δύο ειδών)
Μελετήστε το μοντέλο για τους πληθυσμούς δύο ειδών x, y τα οποία ανταγωνίζονται για το ίδιο είδος τροφής: \begin{align} \frac{d N_1}{dt} & = r_1 N_1\, \left(1 - \frac{N_1}{K_1} - b_{12} \frac{N_2}{K_1}\right) \\ \frac{d N_2}{dt} & = r_2 N_2\, \left(1 - \frac{N_2}{K_2} - b_{21} \frac{N_1}{K_2}\right), \nonumber \end{align} (α) Βρείτε τα σημεία ισορροπίας του συστήματος.
(β) Βρείτε αριθμητική λύση x(t), y(t) (για αρχική συνθήκη της επιλογής σας) και σχεδιάστε τους πληθυσμούς ως συναρτήσεις του χρόνου N₁(t), N₂(t).
(γ) Βρείτε αριθμητικές λύσεις και σχεδιάστε αρκετές καμπύλες στον χώρο των φάσεων (N₁, N₂).

(α) Κάνουμε αδιαστατοποίηση του συστήματος: \begin{align} \frac{d u_1}{dt} & = u_1\, \left(1 - u_1 - a_{12}\,u_2 \right) \\ \frac{d u_2}{dt} & = \rho u_2\, \left(1 - u_2 - a_{21}\,u_1\right), \nonumber \end{align} Έχουμε σημεία ισορροπίας: \[ (u_1, u_2) = (0,0), \quad (u_1, u_2) = (1,0),\quad (u_1, u_2) = (0,1), \] \[ (u_1, u_2) = \left( \frac{1-a_{12}}{1-a_{12} a_{21}}, \frac{1-a_{21}}{1-a_{12} a_{21}} \right). \]

(β) Δείτε έναν κώδικα matlab εδώ και την απαραίτητη function εδώ.

Προσέξτε τα ακόλουθα σημεία:

Επιλέγοντας τις τιμές των παραμέτρων καθορίζουμε και τη θέση του 4ου σημείου ισορροπίας. Επίσης καθορίζει και την ευστάθειά του.
Επιλέγουμε τις αρχικές συνθήκες σε σχέση με τα σημεία ισορροπίας.
Αναρωτηθείτε για την εξέλιξη των $u_1, u_2$ σε σχέση με τα σημεία ισορροπίας.

Ο κώδικας δίνει τις ακόλουθες γραφικές παραστάσεις (για δύο διαφορετικές αρχικές συνθήκες):
two-sppcies solution two-species solution

(γ) Χρειάζεται να κάνουμε μία μικρή αλλαγή στους κώδικες του ερωτήματος (β)...

Ερωτήσεις.

Τι διαδικασία περιγράφει κάθε ένας από τους όρους στα δεξιά μέλη των εξισώσεων του μοντέλου;
Για κάποια μοντέλα έχουμε τον όρο του a₂₁ με το πρόσημο "+". Τι διαδικασίες περιγράφουν αυτά τα μοντέλα;
Βρείτε τα σημεία ισορροπίας και σχεδιάστε το διάγραμμα φάσεων για ένα τέτοιο μοντέλο (με πρόσημο "+" στον όρο του c₂).

Μαθησιακό αποτέλεσμα. Εξοικείωση με λύσεις συστημάτων διαφορικών εξισώσεων και παράστασης λύσεων. Εξοικείωση με ρεαλιστικότερα μοντέλα αλληλεπίδρασης πληθυσμών.

Άσκηση (Μοντέλο επιδημιών SIR)
Μελετήστε το μοντέλο επιδημιών SIR.
(α) Βρείτε τα σημεία ισορροπίας του συστήματος (κάνετε μία δική σας επιλογή για τις σταθερές).
(β) Βρείτε αριθμητική λύση S(t), I(t) (για αρχική συνθήκη της επιλογής σας) και σχεδιάστε τους πληθυσμούς ως συναρτήσεις του χρόνου S(t), I(t).
(γ) Βρείτε αριθμητικές λύσεις και σχεδιάστε αρκετές καμπύλες στον χώρο των φάσεων (S, I).

(α)

(β)

(γ)

Ερωτήσεις.

Τι διαδικασία περιγράφει κάθε ένας από τους όρους στα δεξιά μέλη των εξισώσεων του μοντέλου;
Υπάρχει ένα μοντέλο διάδοσης φημών το οποίο μπορεί να μετασχηματισθεί σε ένα μοντέλο ίδιο με το SIR. Δείτε τι σημαίνουν οι λύσεις του μοντέλου αυτού για την περίπτωση επιδημίας και φήμης αντίστοιχα.

Μαθησιακό αποτέλεσμα. Αριθμητική επίλυση και διάγραμμα φάσεων για μοντέλα με άπειρα σημεία ισορροπίας.

Διαγράμματα φάσεων: γραμμικά συστήματα

Άσκηση (Γραμμικά συστήματα εξισώσεων)
Το γενικό γραμμικό σύστημα δύο εξισώσεων έχει τη μορφή \begin{align} \dot{x} & = a x + b y \\ \dot{y} & = c x + d y. \end{align} Επιλέξτε τιμές για τις σταθερές ώστε το σημείο ισορροπίας να είναι (α) σαγματικό σημείο (β) κόμβος ευσταθής ή ασταθής (γ) σπείρα ευσταθής είτε ασταθής, δεξιόστροφη είτε αριστερόστροφή (δ) κέντρο. Σχεδιάστε τα διαγράμματα φάσεων.
[Υπόδειξη: Δείτε έναν απλό κώδικα matlab και την απαραίτητη function.]

(α) Το σύστημα \begin{align} \dot{x} & = -x - 3 y \\ \dot{y} & = 2 y \end{align} έχει ιδιοτιμές $\lambda_1=-1,\; \lambda_2 = 2$, άρα το σημείο ισορροπίας είναι σαγματικό. Τα ιδιοδιανύσματα του πίνακα του συστήματος είναι $\vec{v}_1=(1,0),\; \vec{v}_2=(1,-1)$.

Σχεδιάστε το διάγραμμα φάσεων.

(β) Το σύστημα \begin{align} \dot{x} & = x - 2 y \\ \dot{y} & = 3 x - 4 y \end{align} έχει ιδιοτιμές $\lambda_1=-1,\; \lambda_2 = -2$, άρα το σημείο ισορροπίας είναι ευσταθής κόμβος. Τα ιδιοδιανύσματα του πίνακα του συστήματος είναι $\vec{v}_1=(1,1),\; \vec{v}_2=(2,3)$.

Το πρόγραμμα matlab για λύση γραμμικού συστήματος βρίσκεται εδώ. Η function για το παραπάνω σύστημα βρίσκεται εδώ.

Ο κώδικας δίνει το ακόλουθο αποτέλεσμα:
node

(γ) Το σύστημα \begin{align} \dot{x} & = -x - 2 y \\ \dot{y} & = 2 x - y. \end{align} έχει ιδιοτιμές $\lambda_{1,2}=-1\pm 2i$, άρα το σημείο ισορροπίας είναι ευσταθής σπείρα.

Σχεδιάστε το διάγραμμα φάσεων.

(δ) Δείτε το σύστημα \begin{align} \dot{x} & = ay \\ \dot{y} & = -x \end{align} όπου a είναι σταθερά.

Ερωτήσεις.

Σχεδιάστε τον ασταθή και ευσταθή υπόχωρο στα διαγράμματα φάσεων για την περίπτωση που έχουμε κόμβο ή σαγματικό σημείο. (Για να το κάνετε αυτό θα χρειαστείτε τα ιδιοδιανύσματα του πίνακα του συστήματος. Αυτά μπορείτε να τα βρείτε με το matlab.)
Γράψτε ένα γραμμικό σύστημα το οποίο παράγει μία δεξιόστροφη σπείρα και ένα άλλο το οποίο παράγει αριστερόστροφη σπείρα. Κάνετε τα αντίστοιχα διαγράμματα φάσεων.
Επαναλάβετε ορισμένα παραδείγματα χρησιμοποιώντας τον κώδικα python για τα διαγράμματα φάσεων.

Μαθησιακό αποτέλεσμα. Διαγράμματα φάσεων για όλες τις βασικές περιπτώσεις σημείων ισορροπίας.

Άσκηση (Ειδικές περιπτώσεις σημείων ισορροπίας)
Σχεδιάστε το διάγραμμα φάσεων για το γραμμικό σύστημα \begin{align} \dot{x} & = x + y \\ \dot{y} & = -x + 3 y. \end{align}

Οι ιδιοτιμές του πίνακα του συστήματος είναι ίσες και είναι λ = 2.

Σχεδιάστε το διάγραμμα φάσεων.

Μαθησιακό αποτέλεσμα. Εξοικείωση με ειδικές περιπτώσεις σημείων ισορροπίας.

Άσκηση (Διάγραμμα φάσεων για κέντρο)
Σχεδιάστε το διάγραμμα φάσεων για το σύστημα \begin{align} \dot{x} & = ay \\ \dot{y} & = -x \end{align} χρησιμοποιώντας το ολοκλήρωμα που μπορείτε να βρείτε. (Χρησιμοποιήστε τιμή για το a της επιλογής σας.)

Σημείο ισορροπίας (x,y) = (0,0).

Υπάρχει ολοκλήρωμα του συστήματος: \[ x^2 + a y^2 = C \] για κάθε $C > 0$.

Σχεδιάζουμε πρώτα τη συνάρτηση (ας επιλέξουμε a=2): \[ f(x,y) = x^2 + a y^2. \] Οι εντολές για το gnuplot βρίσκονται εδώ.

Ακολούθως, σχεδιάζουμε ισοσταθμικές για την παραπάνω συνάρτηση και αυτές μας δίνουν το διάγραμμα φάσεων του συστήματος. Οι εντολές για το gnuplot βρίσκονται εδώ.
Παρατηρήστε ότι C=0 στο σημείο ισορροπίας και C > 0 για κάθε άλλη καμπύλη του διαγράμματος φάσεων.

Ερωτήσεις.

Σχεδιάστε το διάγραμμα φάσεων χρησιμοποιώντας το matlab.
Βρείτε το ολοκλήρωμα που χρησιμοποιήσαμε σε αυτή την άσκηση.
Εξηγήστε γιατί οι ισοσταθμικές καμπύλες της f(x,y) δίνουν τις καμπύλες του διαγράμματος φάσεων.

Μαθησιακό αποτέλεσμα. Σχεδιασμός απλών διαγραμμάτων φάσεων με χρήση πρώτων ολοκληρωμάτων.

Διαγράμματα φάσεων: μη-γραμμικά συστήματα

Άσκηση (απλό εκκρεμές) Σχεδιάστε το διάγραμμα στο χώρο φάσεων για το απλό εκκρεμές:
\[ \frac{d^2\theta}{dt^2} + a \sin(\theta) = 0. \]

Θέτουμε ω=dθ/dt, οπότε έχουμε το εξής σύστημα δύο εξισώσεων
dθ/dt=ω
dω/dt=-a sin(θ)

Σταθερά σημεία: θ=kπ, ω=0

Για να λύσουμε το σύστημα και να σχεδιάσουμε καμπύλες στο χώρο φάσεων χρησιμοποιούμε τον κώδικα matlab εδώ και την αντίστοιχη function εδώ.

Ο κώδικας δίνει το ακόλουθο αποτέλεσμα:
pendulum

2ος τρόπος για να σχεδιάσουμε καμπύλες χώρου φάσεων. Οι καμπύλες του χώρου φάσεων δίνονται από την
(dθ/dt)² - a cos(θ) = E, (E:const.)
Μπορούμε να σχεδιάσουμε τις καμπύλες με εντολές ανάλογες με αυτές στην άσκηση για τον αναρμονικό ταλαντωτή είτε αυτές στην άσκηση για το κέντρο.

Ερωτήσεις.

Κάνετε γραμμικοποίηση γύρω από το σημείο ισορροπίας (π, 0). Δείξτε ότι είναι σαγματικό σημείο και βρείτε τους αναλλοίωτους υποχώρους του γραμμικοποιημένου συστήματος (δηλαδή, τα ιδιοδιανύσματα του γραμμικοποιημένου πίνακα).

Μαθησιακό αποτέλεσμα. Μη-γραμμικά δυναμικά συστήματα στη Μηχανική.

Exercise (Anharmonic oscillator)
(a) Plot the phase portrait for the anharmonic oscillator
d²x/dt² = -x + a x³.
(b) Plot the stable and unstable manifolds for the linearized system around one of the saddle points.

(a) Set y=dx/dt, thus we have the system of equations
dx/dt=y
dy/dt=-x + a x³

Fixed Points:
(i) x=0, y=0
(ii) x=1/sqrt(a), y=0
(iii) x=-1/sqrt(a), y=0

Phase curves are given by
dy/dx = (-x + a x³)/y ==> y²/2 + x²/2 - a x⁴/4 = C

Let us define the function f(x,y)=y²/2 + x²/2 - a x⁴/4
and plot it. Choose, for example a=0.01, so that fixed points are (0,0), (10,0), (-10,0).
The graph of f(x,y) can be plotted using the gnuplot commands:

a=0.1
set xrange[-15:15]
set yrange[-10:10]
splot y**2/2.+x**2/2.-a*x**4/4. w l

Then, make a contour plot of f(x,y) which is equivalent to plotting the phase curves. [explain why.]
In gnuplot we may use this script.
[Note that C=0 at the fixed point (0,0) while the curve with C=25 contains the points (+-10,0).]

Finally, to print your graphs you may use the following commands:

set term post
set output 'anharmonic_phaseportrait.eps'
replot

Here is the result:
Anharmonic oscillator - Phase portrait

(b)

Let us choose a=1. The linearized equations around the fixed point (1,0) are dξ/dt = y, dy/dt = 2ξ, where ξ=x-1. Example for code:

set xrange[0.5:1.5]; set yrange[-0.5:0.5]
set cntrparam levels incremental 0,0.125,0.5
set isosamples 100,100
set title 'Anharmonic oscillator: phase curves (solid), linearized eqns (dashed)' font "Helvetica,16"
splot x**2/2+y**2/2-x**4/4 w l, (x-1)**2/2-y**2/4 w l

Here is the result:
Anharmonic oscillator - Stable and unstable manifolds