ΕΠΙΣΤΗΜΟΝΙΚΟΣ ΥΠΟΛΟΓΙΣΜΟΣ Ι
Ονοματεπώνυμο: Αραβανής Κων/νος Α .Μ.:3628
3η
Εργαστηριακή άσκηση
Σχετικά με το format
της αναφοράς (.htm, .docx) πρέπει να σημειώσω ότι δεν
εμφανίζονται οι εικόνες που περιέχονται σε αυτό με Opera
και Mozilla αλλά μόνο με Internet explorer.
Σε περίπτωση πάντως οποιουδήποτε προβλήματος καλύτερα να γίνει διόρθωση του .docx αρχείου.
0 Χαρακτηριστηκά
Συστήματος
Αρχικά θα
περιγράψουμε το σύστημα στο οποίο τρέξαμε τα πειράματα μας:
Το σύστημα
μας είχε το λειτουργικό Microsoft Windows Vista Business και ο
επεξεργαστής του ήταν ένας Intel Core Duo T2300 στα
1,66GHz.
Είχαμε μία ιεραρχεία μνήμης δύο επιπέδων με το πρώτο επίπεδο να αποτελείτε από
δύο κρυφές μνήμες, η μία για data (L1 DCache) και η
άλλη για inst. (L1 D-Cache) και οι
δύο των 32ΚBytes
και με τα ακόλουθα χαρακτηριστικά: 8 way set-associative, 64-byte line size. Η κύρια
μνήμη (L2 Cache) μας είχε
μέγεθος 2048ΚBytes
και ακριβώς τα ίδια χαρακτηριστικά με τις κρυφές. Όσον αφορά το λογισμικό στο
οποίο εκτελέστηκαν τα προγράμματα μας ήταν η Matlab 7.4
1
Επίλυση διαφορικών εξισώσεων
1.1
Προβλήματα συνοριακών τιμών
1)Στο σημείο
αυτό κλειθήκαμε να υλοποιήσουμε συνάρτηση για πρόβλημα συνοριακώ τιμών που να
βρίσκει το μητρώο συντελεστών Α και το δεξιό μέλος F, ώστε Ax=F.
Ο κώδικας που
χρησιμοποιήσαμε για αυτό είναι ο ακόλουθος:
function [A,F]=f1(b1,
c1, d1, xa, xt, n, eidos, ua, ut) %[A,F]=f1(b1,
c1, d1, xa, xt, n, eidos, ua, ut) %I
sinartisi afti linei provlimata sinoriakon sinthikon %---------------------------------------------------------------- %Dexetai
san orismata tou sintelestes b, c, d pou %edo
emfanizontai me ta onomata b1, c1, d1 kai prokeintai %gia
sinartiseis tou x. Gia to d prepei na isxuei: %d(x)=-u"(x)+b(x)*u'(x)+c(x)*u(x)
me b(x), c(x) %ta b1,
c1 kai u(x) tin sinartissei pou
psaxnoume. Episis san %orismata
dinoume kai to diastima orismou [xa,xt] kai to arithmo %ton
isapexonton komvon n. %An tora
gnorizoume ta: %1)|--Dirichlet-Dirichlet--|
% u(xa),u(xt) thetoume ua=u(xa),ut=u(xt) kai
eidos=1 %2)|--Neumann-Dirichlet--| % u'(xa),u(xt) thetoume ua=u'(xa),ut=u(xt)
kai eidos=2 %3)|--Dirichlet-Neumann--| % u(xa),u'(xt) thetoume ua=u(xa),ut=u'(xt)
kai eidos=3 %4)|--Neumann-Neumann--| % u'(xa),u'(xt) thetoume ua=u'(xa),ut=u'(xt)
kai eidos=4 %5)|--Periodiki
sinartisi sta akra--| % eidos=5, ta ua kai ut edo perisseuoun %---------------------------------------------------------------- %San
apotelesma i sinartisi dinei to mitroo sinteleston A kai to %deksio
melos F oste A*x=F %---------------------------------------------------------------- %euresi
ton sinteleston os sinartiseis tou x b=inline(b1,'x'); c=inline(c1,'x'); d=inline(d1,'x'); %euresi
tou diastimatos h pou apexei o enas komvos apo ton allo h=(xt-xa)/(n+1); %parakato
opou emfanizontai ta aa, bb, cc ennoountai ta α, β, γ %apo tin
theoria %elegxos
gia to se poia periptosi vriskomaste %Dirichlet-Dirichlet if eidos==1 aa=zeros(n,1); bb=zeros(n-1,1); cc=zeros(n-1,1); F=zeros(n,1); %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] for i=1:1:n-1 aa(i,1)=2/(h^2) + c(xa+i*h); bb(i,1)=-1/(h^2) + b(xa+i*h)/(2*h); cc(i,1)=-1/(h^2) -b(xa+(i+1)*h)/(2*h); end aa(n,1)=2/(h^2) + c(xa+n*h); %euresi tou F F(1,1)=d(xa+h)+ua/(h^2)+b(xa+h)*ua/(2*h); for i=2:1:n-1 F(i,1)=d(xa+i*h); end
F(n,1)=d(xa+n*h)+ut/(h^2)-b(xa+n*h)*ut/(2*h); end %Neumann-Dirichlet if eidos==2 aa=zeros(n+1,1); bb=zeros(n,1); cc=zeros(n,1); F=zeros(n+1,1); %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] aa(1,1)=2/(h^2)+c(xa); bb(1,1)=-2/(h^2); cc(1,1)=-1/(h^2) -b(xa+h)/(2*h); for i=2:1:n aa(i,1)=2/(h^2)+c(xa+(i-1)*h); bb(i,1)=-1/(h^2) +
b(xa+(i-1)*h)/(2*h); cc(i,1)=-1/(h^2) -b(xa+i*h)/(2*h); end
aa(n+1,1)=2/(h^2)+c(xa+n*h); %euresi tou F F(1,1)=d(xa)-b(xa)*ua-2*ua/h; for i=2:1:n F(i,1)=d(xa+(i-1)*h); end
F(n+1,1)=d(xa+n*h)+ut/(h^2)-b(xa+n*h)*ut/(2*h); end %Dirichlet-Neumann if eidos==3 aa=zeros(n+1,1); bb=zeros(n,1); cc=zeros(n,1); F=zeros(n+1,1); %euresi ton α, β, γ pou xreiazontai
gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] for i=1:1:n-1 aa(i,1)=2/(h^2)+c(xa+i*h); bb(i,1)=-1/(h^2) + b(xa+i*h)/(2*h); cc(i,1)=-1/(h^2)
-b(xa+(i+1)*h)/(2*h); end aa(n,1)=2/(h^2)+c(xa+n*h); bb(n,1)=-1/(h^2) + b(xa+n*h)/(2*h); cc(n,1)=-2/(h^2); aa(n+1,1)=2/(h^2)+c(xa+(n+1)*h); %euresi tou F F(1,1)=d(xa+h)+ua/(h^2)+b(xa+h)*ua/(2*h); for i=2:1:n F(i,1)=d(xa+i*h); end F(n+1,1)=d(xt)-b(xt)*ut+2*ut/h; end %Neumann-Neumann if eidos==4 aa=zeros(n+2,1); bb=zeros(n+1,1); cc=zeros(n+1,1); F=zeros(n+2,1); %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] aa(1,1)=2/(h^2)+c(xa); bb(1,1)=-2/(h^2); cc(1,1)=-1/(h^2) -b(xa+h)/(2*h); for i=2:1:n aa(i,1)=2/(h^2)+c(xa+(i-1)*h); bb(i,1)=-1/(h^2) +
b(xa+(i-1)*h)/(2*h); cc(i,1)=-1/(h^2) -b(xa+i*h)/(2*h); end aa(n+1,1)=2/(h^2)+c(xa+n*h); bb(n+1,1)=-1/(h^2) + b(xa+n*h)/(2*h); cc(n+1,1)=-2/(h^2); aa(n+2,1)=2/(h^2)+c(xt); %euresi tou F F(1,1)=d(xa)-b(xa)*ua-2*ua/h; for i=2:1:n+1 F(i,1)=d(xa+(i-1)*h); end F(n+2,1)=d(xt)-b(xt)*ut+2*ut/h; end %Periodiki
sinartisi sta akra if eidos==5 aa=zeros(n+1,1); bb=zeros(n,1); cc=zeros(n,1); F=zeros(n+1,1); %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] for i=1:1:n aa(i,1)=2/(h^2) + c(xa+i*h); bb(i,1)=-1/(h^2) + b(xa+i*h)/(2*h); cc(i,1)=-1/(h^2)
-b(xa+(i+1)*h)/(2*h); end aa(n+1,1)=2/(h^2) + c(xt); %euresi tou F for i=1:1:n+1 F(i,1)=d(xa+i*h); end end %dimiourgia
tou tridiagoniou mitroou A=trid[γ,β,α] A=GALLERY('tridiag',cc,aa,bb); %Periodiki
sinartisi sta akra if eidos==5 %se afti ti peritosei prepei na
prosthesoume alla dio stoixeia %ston mexri stigmis tridiagonio
A pou dilonoun pera ton allon %oti ta u(xa) kai u(xt) einai
isa A(1,n+1)=-1/(h^2)-b(xa+h)/(2*h); A(n+1,1)=-1/(h^2)+b(xt)/(2*h); end |
Στο σημείο αυτό
ελένξαμε τον αλγοριθμό μας για την ορθότητα του για κάθε μια από τις
περιπτώσεις ξεχωριστά:
Α) Dirichlet Dirichlet à eidos=1
·
Ψάχνουμε
για την συνάρτηση u(x)=x^3
·
b(x)=sin(x)
·
c(x)=x+2
· d(x)=-u(x)+b(x)*u(x)+c(x)*u(x)
· xa=-10 à u(xa)=-1000
· xt=10 à u(xt)=1000
· n=10000
>> [A,F]=f1('sin(x)','x+2','-6*x+sin(x)*(3*(x^2))+(x+2)*(x^3)',-10,10,10000,1,-1000,1000); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
B) Neumann Dirichlet à eidos=2
·
Ψάχνουμε
για την συνάρτηση u(x)=x^3
·
b(x)=sin(x)
·
c(x)=x+2
· d(x)=-u(x)+b(x)*u(x)+c(x)*u(x)
· xa=-10 à u(xa)=300
· xt=10 à u(xt)=1000
· n=10000
>> [A,F]=f1('sin(x)','x+2','-6*x+sin(x)*(3*(x^2))+(x+2)*(x^3)',-10,10,10000,2,300,1000); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
Γ) Dirichlet Neumannà eidos=3
·
Ψάχνουμε
για την συνάρτηση u(x)=x^3
·
b(x)=sin(x)
·
c(x)=x+2
· d(x)=-u(x)+b(x)*u(x)+c(x)*u(x)
· xa=-10 à u(xa)=-1000
· xt=10 à u(xt)=300
· n=10000
>> [A,F]=f1('sin(x)','x+2','-6*x+sin(x)*(3*(x^2))+(x+2)*(x^3)',-10,10,10000,3,-1000,300); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
Δ) Neumann Neumannà eidos=4
·
Ψάχνουμε
για την συνάρτηση u(x)=x^3
·
b(x)=sin(x)
·
c(x)=x+2
· d(x)=-u(x)+b(x)*u(x)+c(x)*u(x)
· xa=-10 à u(xa)=300
· xt=10 à u(xt)=300
· n=10000
>> [A,F]=f1('sin(x)','x+2','-6*x+sin(x)*(3*(x^2))+(x+2)*(x^3)',-10,10,10000,4,300,300); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
Ε)Περιοδική
συνάρτηση με ακρα δοσμένα τα xa, xa+T à eidos=5
·
Ψάχνουμε
για την συνάρτηση u(x)=x^3
·
b(x)=sin(x)
·
c(x)=x+2
· d(x)=-u(x)+b(x)*u(x)+c(x)*u(x)
· xa=0
· xt=2π
· n=10000
>> [A,F]=f1('cos(x)','x','sin(x)+cos(x)*(cos(x))+(x)*(sin(x))', 0, 2*pi,10000, 5); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
2,3) Στο
σημείο αυτό της άσκησης θα αποδείξουμε θεωρητικά πότε ένα μητρώο που εχει
προκύψει από συνοριακές τιμές
Dirichlet- Dirichlet, Neumann-Dirichlet, Dirichlet- Neumann, Neumann- Neumann ή Περιοδικές είναι συμμετρικό, θετικά
ορισμένο ή ημιορισμένο, Toeplitz, αραιό, κυκλοτερές και πότε μη αντιστρέψιμο.
Ένα μητρώο Dirichlet- Dirichlet, Neumann-Dirichlet, Dirichlet- Neumann ή Neumann- Neumann έχει τη μορφή
Α=|α1 β1 0 0 ... 0 0 0 |
|γ2 α2 β2 0 ...
0 0 0
|
|0
γ3 α3 β3 ... 0 0
0 |
|..........................................|
|0
0 0 0 ...
γκ-1 ακ-1 βκ-1 |
|0
0 0 0 ...
0 γκ ακ |
(Για Dirichlet- Dirichlet κ=n, για Neumann-Dirichlet και Dirichlet- Neumann κ=n+1 και για Neumann- Neumann κ=n+2)
Α) Για να
είναι συμμετρικό ένα τέτοιο μητρώο θα πρέπει να ισχύει γj+1=βj =>
-1/(h^2) - bj+1/(2*h)= -1/(h^2) + bj/(2*h) =>
- bj+1 = bj =>
bj=0 |
Β) Για να
είναι θετικά ορισμένο ή ημιορισμένο ένα τέτοιο μητρώο θα πρέπει να ισχύει: xAx>0 για κάθε μη μηδενικό διάνυσμα x.
Γ) Για να είναι Toeplitz ένα τέτοιο μητρώο θα πρέπει να ισχύει ότι:
αj=ct, βj=ct και γj=ct
βj= βj+1
γj= γj+1
-1/(h^2) + bj/(2*h) =-1/(h^2) + bj+1/(2*h)
-1/(h^2) - bj/(2*h) =-1/(h^2) bj+1/(2*h)
bj = bj+1
Δ) Για να
είναι Αραιό ένα τέτοι μητρώο θα πρέπει να ισχύει ότι μόνο το 20% των στοιχειων
του είναι μη μηδενικά το πολύ.
Δηλαδή στη
περίπτωση μας που έχουμε 3κ (~=κ+(κ-1)+(κ-1)) μη μηδενικά στοιχεία θα πρέπει να
ισχύει:
3κ/κ^2 <
20% =>
κ> 15 |
Δηλαδή για:
·
Dirichlet Dirichlet
n>15 |
·
Neumann Dirichlet και Dirichlet - Neumann
n>14 |
·
Neumann - Neumann
n>13 |
E) Είναι
προφανές ότι για ένα τέτοιο μητρώο δεν μπορεί να είναι κυκλοτρεπές παρά μόνον
αν:
γj=0
-1/(h^2) + bj/(2*h) =0
-1/(h^2) - bj/(2*h) =0
Πράγμα άτοπο
άρα δεν μπορεί να είναι κυκλοτρεπές σε καμμιά περίπτωση
Ζ) Για να
είναι μη αντιστρέψιμο ένα τέτοιο μητρώο πρέπει να υπάρχει τουλάχιστον μία
ιδιοτιμή του που να είναι 0.
Ένα μητρώο
περιοδικών έχει τη μορφή:
Α=|α1 β1 0 0 ...
0 0 γ1 |
|γ2 α2 β2 0 ...
0 0 0
|
|0
γ3 α3 β3 ... 0 0
0 |
|...........................................|
|0
0 0 0 ...
γn-1
αn-1 βn-1 |
|βn
0 0 0 ...
0 γn αn |
Α) Για να
είναι συμμετρικό ένα τέτοιο μητρώο θα πρέπει να ισχύει γj+1=βj και γ1=βn
Μόνο για το πώτο
σκέλος δηλαδή θα πρέπει να ισχύει
-1/(h^2) - bj+1/(2*h)= -1/(h^2) + bj/(2*h) =>
- bj+1 = bj =>
bj=0
που
συνεπάγεται άμεσα και το δεύτερο σκέλος
Επομένος ένα
τέτοιο μητρώο είναι συμμετρικό εάν:
bj=0 |
Β) Για να
είναι θετικά ορισμένο ή ημιορισμένο ένα τέτοιο μητρώο θα πρέπει να ισχύει: xAx>0 για κάθε μη μηδενικό διάνυσμα x.
Γ) Για να είναι Toeplitz ένα τέτοιο μητρώο θα πρέπει να ισχύει ότι:
αj=ct, βj=ct και γj=ct
βj= βj+1
γj= γj+1
-1/(h^2) + bj/(2*h) =-1/(h^2) + bj+1/(2*h)
-1/(h^2) - bj/(2*h) =-1/(h^2) bj+1/(2*h)
bj = bj+1
Δ) Για να
είναι Αραιό ένα τέτοι μητρώο θα πρέπει να ισχύει ότι μόνο το 20% των στοιχειων
του είναι μη μηδενικά το πολύ.
Δηλαδή στη
περίπτωση μας που έχουμε 3n
(=n+(n-1) + (n-1) + 1 + 1) μη μηδενικά στοιχεία θα πρέπει να ισχύει:
3n/n^2 < 20% =>
n> 15 |
E) Είναι
προφανές ότι για ένα τέτοιο μητρώο δεν μπορεί να είναι κυκλοτρεπές παρά μόνον
αν:
βj= βj+1
γj= γj+1
-1/(h^2) + bj/(2*h) =-1/(h^2) + bj+1/(2*h)
-1/(h^2) - bj/(2*h) =-1/(h^2) bj+1/(2*h)
bj = bj+1
Ζ) Για να
είναι μη αντιστρέψιμο ένα τέτοιο μητρώο πρέπει να υπάρχει τουλάχιστον μία
ιδιοτιμή του που να είναι 0.
Ένας
συνοπτικός πίνακας που περιγράφει με ποιές προϋποθέσης το Α μπορεί να είναι
συμμετρικό, θετικά ορισμένο ή ημιορισμένο, Toeplitz, αραιό, κυκλοτερές ή μη αντιστρέψιμο φαίνεται
παρακάτω:
|
Dirichlet- Dirichlet |
Neumann-Dirichlet |
Dirichlet- Neumann |
Neumann- Neumann |
Περιοδικές |
Συμμετρικό |
bj=0 |
bj=0 |
bj=0 |
bj=0 |
bj=0 |
Θετικά
ορισμένο ή ημιορισμένο |
xAx>0 για κάθε μη μηδενικό διάνυσμα x |
xAx>0 για κάθε μη μηδενικό διάνυσμα x |
xAx>0 για κάθε μη μηδενικό διάνυσμα x |
xAx>0 για κάθε μη μηδενικό διάνυσμα x |
xAx>0 για κάθε μη μηδενικό διάνυσμα x |
Toeplitz |
cj=cj+1 bj =bj+1 |
cj=cj+1 bj =bj+1 |
cj=cj+1 bj =bj+1 |
cj=cj+1 bj =bj+1 |
cj=cj+1 bj =bj+1 |
Αραιό |
n>15 |
n>14 |
n>14 |
n>13 |
n>15 |
Κυκλοτερές |
Δεν μπορεί
να είναι |
Δεν μπορεί να είναι |
Δεν μπορεί να είναι |
Δεν μπορεί να είναι |
cj=cj+1 bj =bj+1 |
Μη αντιστρέψιμο |
Τουλάχιστον
μια ιδιοτιμή του Α να είναι ίση με 0 |
Τουλάχιστον
μια ιδιοτιμή του Α να είναι ίση με 0 |
Τουλάχιστον
μια ιδιοτιμή του Α να είναι ίση με 0 |
Τουλάχιστον
μια ιδιοτιμή του Α να είναι ίση με 0 |
Τουλάχιστον
μια ιδιοτιμή του Α να είναι ίση με 0 |
Για να ελένξουμε
τα παραπάνω ο κώδικας που φτιάξαμε είναι ο ακόλουθος:
function
[A,flags]=f2(b1, c1, d1, xa, xt, n, eidos, ua, ut) %[flags]=f2(b1,
c1, d1, xa, xt, n, eidos, ua, ut) %I
sinartisi afti vriskei to eidos tou mitroou pou prokiptei apo %ena
provlima sinoriakon sinthikon %---------------------------------------------------------------- %Dexetai
san orismata tou sintelestes b, c, d pou %edo
emfanizontai me ta onomata b1, c1, d1 kai prokeintai %gia
sinartiseis tou x. Gia to d prepei na isxuei: %d(x)=-u"(x)+b(x)*u'(x)+c(x)*u(x)
me b(x), c(x) %ta b1,
c1 kai u(x) tin sinartissei pou
psaxnoume. Episis san %orismata
dinoume kai to diastima orismou [xa,xt] kai to arithmo %ton
isapexonton komvon n. %An tora
gnorizoume ta: %1)|--Dirichlet-Dirichlet--|
% u(xa),u(xt) thetoume ua=u(xa),ut=u(xt) kai
eidos=1 %2)|--Neumann-Dirichlet--| % u'(xa),u(xt) thetoume ua=u'(xa),ut=u(xt)
kai eidos=2 %3)|--Dirichlet-Neumann--| % u(xa),u'(xt) thetoume ua=u(xa),ut=u'(xt)
kai eidos=3 %4)|--Neumann-Neumann--| % u'(xa),u'(xt) thetoume ua=u'(xa),ut=u'(xt)
kai eidos=4 %5)|--Periodiki
sinartisi sta akra--| % eidos=5, ta ua kai ut edo perisseuoun %---------------------------------------------------------------- %San
apotelesma i sinartisi dinei ena pinaka flags pou ta stoixeia tou mas %deixnoun
an to prokipton mitroo einai %1)simmetriko--->flags(1,1)=1 %2)thetika
orismeno i imiorismeno--->flags(2,1)=1 %3)toeplitz--->flags(3,1)=1 %4)araio--->flags(4,1)=1 %5)kikloteres--->flags(5,1)=1 %6)mi
antistrepsimo--->flags(6,1)=1 %---------------------------------------------------------------- %euresi
ton sinteleston os sinartiseis tou x b=inline(b1,'x'); c=inline(c1,'x'); d=inline(d1,'x'); %euresi
tou diastimatos h pou apexei o enas komvos apo ton allo h=(xt-xa)/(n+1); %parakato
opou emfanizontai ta aa, bb, cc ennoountai ta α, β, γ %apo tin
theoria %elegxos
gia to se poia periptosi vriskomaste %Dirichlet-Dirichlet if eidos==1 aa=zeros(n,1); bb=zeros(n-1,1); cc=zeros(n-1,1); k=n; %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] for i=1:1:n-1 aa(i,1)=2/(h^2) + c(xa+i*h); bb(i,1)=-1/(h^2) + b(xa+i*h)/(2*h); cc(i,1)=-1/(h^2)
-b(xa+(i+1)*h)/(2*h); end aa(n,1)=2/(h^2) + c(xa+n*h); end %Neumann-Dirichlet if eidos==2 aa=zeros(n+1,1); bb=zeros(n,1); cc=zeros(n,1); k=n+1; %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] aa(1,1)=2/(h^2)+c(xa); bb(1,1)=-2/(h^2); cc(1,1)=-1/(h^2) -b(xa+h)/(2*h); for i=2:1:n aa(i,1)=2/(h^2)+c(xa+(i-1)*h); bb(i,1)=-1/(h^2) +
b(xa+(i-1)*h)/(2*h); cc(i,1)=-1/(h^2) -b(xa+i*h)/(2*h); end
aa(n+1,1)=2/(h^2)+c(xa+n*h); end %Dirichlet-Neumann if eidos==3 aa=zeros(n+1,1); bb=zeros(n,1); cc=zeros(n,1); k=n+1; %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] for i=1:1:n-1 aa(i,1)=2/(h^2)+c(xa+i*h); bb(i,1)=-1/(h^2) + b(xa+i*h)/(2*h); cc(i,1)=-1/(h^2)
-b(xa+(i+1)*h)/(2*h); end aa(n,1)=2/(h^2)+c(xa+n*h); bb(n,1)=-1/(h^2) + b(xa+n*h)/(2*h); cc(n,1)=-2/(h^2); aa(n+1,1)=2/(h^2)+c(xa+(n+1)*h); end %Neumann-Neumann if eidos==4 aa=zeros(n+2,1); bb=zeros(n+1,1); cc=zeros(n+1,1); k=n+2; %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] aa(1,1)=2/(h^2)+c(xa); bb(1,1)=-2/(h^2); cc(1,1)=-1/(h^2) -b(xa+h)/(2*h); for i=2:1:n aa(i,1)=2/(h^2)+c(xa+(i-1)*h); bb(i,1)=-1/(h^2) +
b(xa+(i-1)*h)/(2*h); cc(i,1)=-1/(h^2) -b(xa+i*h)/(2*h); end aa(n+1,1)=2/(h^2)+c(xa+n*h); bb(n+1,1)=-1/(h^2) + b(xa+n*h)/(2*h); cc(n+1,1)=-2/(h^2); aa(n+2,1)=2/(h^2)+c(xt); end %Periodiki
sinartisi sta akra if eidos==5 aa=zeros(n+1,1); bb=zeros(n,1); cc=zeros(n,1); k=n+1; %euresi ton α, β, γ pou
xreiazontai gia tin dimiourgia tou %tridiagoniou mitroou
A=trid[γ,β,α] for i=1:1:n aa(i,1)=2/(h^2) + c(xa+i*h); bb(i,1)=-1/(h^2) + b(xa+i*h)/(2*h); cc(i,1)=-1/(h^2)
-b(xa+(i+1)*h)/(2*h); end aa(n+1,1)=2/(h^2) + c(xt); end %dimiourgia
tou tridiagoniou mitroou A=trid[γ,β,α] A=GALLERY('tridiag',cc,aa,bb); %Periodiki
sinartisi sta akra if eidos==5 %se afti ti peritosei prepei na
prosthesoume alla dio stoixeia %ston mexri stigmis tridiagonio
A pou dilonoun pera ton allon %oti ta u(xa) kai u(xt) einai
isa A(1,n+1)=-1/(h^2)-b(xa+h)/(2*h); A(n+1,1)=-1/(h^2)+b(xt)/(2*h); end flags=zeros(6,1); %an to
b(x)=0 tote to A einai simmetriko if b1=='0' flags(1,1)=1; end x=rand(k,1); %an to
x'*A*x>0 tote to A einai thetika orismeno i imiorismeno if
((x')*A*x)>0 flags(2,1)=1; end notequal=1; for i=1:1:k-2 %ta antidiametrika stoixeia tis
diagoniou einai isa? if A(i+1,i)~=A(i+2,i+1) notequal=0; end if A(i,i+1)~=A(i+1,i+2) notequal=0; end
if A(i,i)~=A(i+1,i+1) notequal=0; end end if
A(k-1,k-1)~=A(k,k) notequal=0; end %einai
toeplitz? if notequal==1 flags(3,1)=1; %einai kikloteres? if eidos==5 flags(5,1)=1; end end %an
k>15 einai araio if k>15 flags(4,1)=1; end %euresi
ton idiotimon tou A [V]=eigs(A); zeroidiotimi=0; for
i=1:1:length(V) %einai isi me 0 i idiotimi se
afti ti thesi tou dianismatos? if V(i,1)==0 zeroidiotimi=1; end end %an
iparxei mideniki idiotimi tote to A einai mi antistrapsimo if
zeroidiotimi==1 flag(6,1)=1; end |
Αν για
παράδειγμα δώσουμε σαν είσοδο στη συνάρτηση αυτή το ακόλουθο:
>> [A,flags]=f2('sin(x)','x+2','-6*x+sin(x)*(3*(x^2))+(x+2)*(x^3)',-10,10,500,1,-1000,1000) |
Παίρνουμε σαν
απάντηση το ακόλουθο:
flags = 0 1 0 1 0 0 |
Που σημαίνει
ότι το Α
·
Δεν
είναι συμμετρικό
·
Είναι
θετικά ορισμένο ή ημιορισμένο
·
Είναι
Toeplitz
·
Είναι
αραιό
·
Δεν
είναι κυκλοτερές
·
Δεν
είναι μη αντιστρέψιμο
Ένώ αν δώσουμε
το ακόλουθο:
>> [A,flags]=f2('sin(x)','x+2','-6*x+sin(x)*(3*(x^2))+(x+2)*(x^3)',-10,10,500,5) |
Παίρνουμε σαν
απάντηση το ακόλουθο:
flags = 0 1 0 1 1 0 |
Που σημαίνει
ότι το Α
·
Είναι
συμμετρικό
·
Είναι
θετικά ορισμένο ή ημιορισμένο
·
Είναι
Toeplitz
·
Είναι
αραιό
·
Eίναι
κυκλοτερές
·
Δεν
είναι μη αντιστρέψιμο
Ένα δένδρο απόφασης που μπορεί να συνοψίσει όλες τις παραπάνω δυνατές
περιπτώσεις στις οποίες μπορεί να
βρισκεται ένα μητρώο Α σε προβλήματα συνοριακών τιμών οπου για τα ακρα μπορεί
να έχουμε συνοριακές τιμές Dirichlet- Dirichlet, Neumann-Dirichlet, Dirichlet- Neumann, Neumann- Neumann ή Περιοδικές
φαίνεται παρακάτω (Η ιδέα αυτή για την αναπαράσταση του συνδιασμού τον
χαρακτηριστικών του μητρώου Α πάρθηκε από το http://scgroup.hpclab.ceid.upatras.gr/class/SC/Lectures/Set17.pdf σημειώσεις του κ.Γαλλόπουλου (σελ.11))
Για
την εκτύπωση των γραφημάτων φτιάξαμε την ακόλουθη συνάρτηση που θεωρεί σαν
δεδομένο ότι η συνάρτηση που ψάχνουμε είναι η
sin(x) :
function f3(j) %f3(j),
j=1:1:5 %i
sinartisi afti dexetai san orismata akeraious apo to 1 eos to 5 oste %o
xristis na epileksi metaksi %1)|--Dirichlet-Dirichlet--|
%2)|--Neumann-Dirichlet--| %3)|--Dirichlet-Neumann--| %4)|--Neumann-Neumann--| %5)|--Periodiki
sinartisi sta akra--| %kai
vriskei to deikti katastasis tou A gia n=2^i me i=6:1:10 %to A
mporei na exei diafora xaraktiristika %i
sinartisi mas afti ftiaxnei 4on eidon A me diaforetiko sindiasmo %xaraktiristikon,
ta eidoi ton A einai ta A,B,C,D pou fainontai sta sxolia %esoterika
tis sinartisis gia to poia ta xaraktiristika tou kathe mitroou %gia
kathe provlima sinoriakon timon condition=zeros(4,5); n=[2^6;2^7;;2^8;2^9;2^10]; d=[1 0 0;2 1
0;3 0 1;4 1 1;5 0 0]; %A %thetika
orismeno-araio---->Dirichlet-Dirichlet %thetika
orismeno-araio---->Neumann-Dirichlet %thetika
orismeno-araio---->Dirichlet-Neumann %thetika
orismeno-araio---->Neumann-Neumann %thetika
orismeno-araio---->Periodikes for i=6:1:10 [A,F]=f1('cos(x)','x','sin(x)+cos(x)*(cos(x))+(x)*(sin(x))', 0, 2*pi,2^i,
d(j,1),d(j,2),d(j,3)); condition(1,i-5)=condest(A); end %B %simmetriko-thetika
orismeno-araio---->Dirichlet-Dirichlet %simmetriko-thetika
orismeno-araio---->Neumann-Dirichlet %simmetriko-thetika
orismeno-araio---->Dirichlet-Neumann %simmetriko-thetika
orismeno-araio---->Neumann-Neumann %simmetriko-thetika
orismeno-araio---->Periodikes for i=6:1:10 [A,F]=f1('0','x','sin(x)+0*(cos(x))+(x)*(sin(x))', 0, 2*pi,2^i,
d(j,1),d(j,2),d(j,3)); condition(2,i-5)=condest(A); end %C %thetika
orismeno-toeplitz-araio---->Dirichlet-Dirichlet %thetika
orismeno-toeplitz-araio---->Neumann-Dirichlet %thetika
orismeno-toeplitz-araio---->Dirichlet-Neumann %thetika
orismeno-toeplitz-araio---->Neumann-Neumann %thetika
orismeno-toeplitz-araio-kikloteres---->Periodikes for i=6:1:10
[A,F]=f1('5','1','sin(x)+5*(cos(x))+(1)*(sin(x))', 0, 2*pi,2^i,
d(j,1),d(j,2),d(j,3)); condition(3,i-5)=condest(A); end %D %simmetriko-thetika
orismeno-toeplitz-araio---->Dirichlet-Dirichlet %simmetriko-thetika
orismeno-toeplitz-araio---->Neumann-Dirichlet %simmetriko-thetika
orismeno-toeplitz-araio---->Dirichlet-Neumann %simmetriko-thetika
orismeno-toeplitz-araio---->Neumann-Neumann %simmetriko-thetika
orismeno-toeplitz-araio-kikloteres---->Periodikes for i=6:1:10 [A,F]=f1('0','1','sin(x)+0*(cos(x))+(1)*(sin(x))', 0, 2*pi,2^i,
d(j,1),d(j,2),d(j,3)); condition(4,i-5)=condest(A); end semilogy(n,
condition(1,:), 'r.-'); %
kanoume hold gia na emfanistoun se koini grafiki parastasi hold on; semilogy(n,
condition(2,:), 'm.-'); semilogy(n,
condition(3,:), '.-'); semilogy(n,
condition(4,:), 'g.-'); legend('A','B','C','D') hold; |
Τα γραφήματα
που παίρνουμε σαν αποτέλεσμα θέτοντας j=1:1:5 είναι τα ακόλουθα:
A)Βάζουμε άκρα Dirichlet-Dirichlet
>> f3(1) |
Το Α μητρωο
που αναπαριστάτε από κάθε γραμμή έχει τα ακόλουθα χαρακτηριστηκά
Α: Θετικά ορισμενο/Αραιό
Β: Συμμετρικό/ Θετικά ορισμενο/Αραιό
C: Θετικά ορισμενο/Toeplitz/Αραιό
D: Συμμετρικό /Θετικά ορισμενο/Toeplitz/Αραιό
Β) Βάζουμε άκρα
Neumann-Dirichlet
>> f3(2) |
Το Α μητρωο
που αναπαριστάτε από κάθε γραμμή έχει τα ακόλουθα χαρακτηριστηκά
Α: Θετικά ορισμενο/Αραιό
Β: Συμμετρικό/ Θετικά ορισμενο/Αραιό
C: Θετικά ορισμενο/Toeplitz/Αραιό
D: Συμμετρικό /Θετικά ορισμενο/Toeplitz/Αραιό
Γ) Βάζουμε άκρα
Dirichlet-Neumann
>> f3(3) |
Το Α μητρωο
που αναπαριστάτε από κάθε γραμμή έχει τα ακόλουθα χαρακτηριστηκά
Α: Θετικά ορισμενο/Αραιό
Β: Συμμετρικό/ Θετικά ορισμενο/Αραιό
C: Θετικά ορισμενο/Toeplitz/Αραιό
D: Συμμετρικό /Θετικά ορισμενο/Toeplitz/Αραιό
Δ) Βάζουμε άκρα Neumann- Neumann
>> f3(4) |
Το Α μητρωο
που αναπαριστάτε από κάθε γραμμή έχει τα ακόλουθα χαρακτηριστηκά
Α: Θετικά ορισμενο/Αραιό
Β: Συμμετρικό/ Θετικά ορισμενο/Αραιό
C: Θετικά ορισμενο/Toeplitz/Αραιό
D: Συμμετρικό /Θετικά ορισμενο/Toeplitz/Αραιό
Ε) Βάζουμε άκρα περιοδικά
>> f3(5) |
Επειδή δεν
φαινεται καλά εδώ για το που είναι η C κάναμε zoom στη Matlab και φάνηκε ότι έχει συμπέσει με την D.
Το Α μητρωο
που αναπαριστάτε από κάθε γραμμή έχει τα ακόλουθα χαρακτηριστηκά
Α: Θετικά ορισμενο/Αραιό
Β: Συμμετρικό/ Θετικά ορισμενο/Αραιό
C: Θετικά ορισμενο/Toeplitz/Αραιό/ Κυκλοτερές
D: Συμμετρικό /Θετικά ορισμενο/Toeplitz/Αραιό/Κυκλοτερές
Από τα
παραπάνω είναι φανερό ότι όσο αυξάνεται το μέγεθος του μητρώου Α (λόγω του ότι
αυξάνεται το n) αυξάνεται και
ο δείκτης κατάστασης. Αυτό ισχύει οπως
γίνεται σαφές από όλα τα παραπάνω σχήματα και για Dirichlet- Dirichlet, Neumann-Dirichlet, Dirichlet- Neumann, Neumann- Neumann και για
Περιοδικά.
5) Στο σημείο
αυτό εξετάζω τη λύση της διαφορικής εξίσωσης ως προς τη μονοτονία για τις δύο
παρακάτω περιπτώσεις:
Α) Για c(x)=d(x)=0, στο διάστημα [0,1] με συνοριακές
τιμές u(0)=0, u(1)=1, b(x)=100 και μεταβλητό h=1/(κ+1) ή με άλλα λόγια μεταβλητό κ
και
Β) Για c(x)=d(x)=0, στο διάστημα [0,1] με συνοριακές
τιμές u(0)=0, u(1)=1, b(x)=ΑΜ-ΕΓ=3628-1987=1641
(ΑΜ=Αριθμός Μητρώου, ΕΓ=Έτος Γέννησης) και μεταβλητό h=1/(κ+1) ή με άλλα λόγια μεταβλητό κ
Παρακάτω
φαίνονται τα αποτελέσματα που λάβαμε για την αναζητούμενη διαφορική εξίσωση σε
κάθε μία από τις 2 περιπτώσεις. Έδώ πρέπει να σημειώσουμε ότι οι διαφορικές
εξισώσεις που αναζητήσαμε για την απόδειξη των σκέψεων που θα καταθέσουμε
παρακάτω ήταν πολύ περισσότερες ώστε να είμαστε σίγουροι για την ορθότητα τους.
Η κατακλείδα
μας σχετικά με το ερώτημα αυτού του μέρους της άσκησης είναι ότι όταν το κ που
θέτουμε στη συναρτηση είναι μικροτερο από
το b(x)/2 δέν μπορούμε να αποφανθούμε για την
μονοτονία της διαφορικής εξίσωσης καθώς αυτή μεταβάλλεται καθώς μεταβάλλεται
και το κ (Συγκεκριμενα όσο πιο μικρό το κ τοσο πιο πολλές είναι οι αλλαγές στη
μονοτονία της συναρτησης στα διάφορα διαστήματα). Επίσης όταν θέτουμε στην f1 κ το οποίο να είναι μεγαλύτερο του b(x)/2 η διαφορική εξίσωση που παίρνουμε σαν απάντηση είναι αυξουσα (για οποιεσδήποτε
τιμές του κ στο συγκεκριμένο διάστημα).Τέλος όταν το κ παίρνει την τιμή b(x)/2 παρατηρούμε στην πρώτη περιπτωση τη διαφορική που προκύπτει να είναι
παράλληλη στον άξονα x και
στη δεύτερη περίπτωση παράλληλη στον y,δηλαδή θα μπορούσε να χαρακτηριστεί ίσως και πάλι αυξουσα.
Τα
συμπεράσματά μας μπορούν να συνοψιστούν στον παρακάτω πίνακα:
κ<b(x)/2 |
κ>=b(x)/2 |
Το πεδίο
ορισμού της διαφορικής εξίσωσής μας
σπάει σε περισσότερα από ένα διαστήματα στα οποία η μονοτονία της διαφορικής
εξίσωσης αλλαζει απο αυξουσα σε φθίνουσα και ανάποδα. Όσο
μικρότερο είναι το κ τόσο πιο έντονο είναι αυτό το φαινόμενο |
Η διαφορική
εξίσωση είναι αύξουσα. |
Μερικές από τα
συμαντικότερα αποτελέσματα που λάβαμε από την f1 για κάθε περίπτωση φαίνονται παρακάτω.
Συγκεκριμάνα παρουσιάζουμε τις εξής περιπτώσεις
·
κ<< b(x)/2
·
κ<b(x)/2
·
κ~= b(x)/2
·
κ> b(x)/2
·
κ>> b(x)/2
Α)
Ι)κ<< b(x)/2
>> [A,F]=f1('100','0','0', 0, 1, 5, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
ΙΙ)κ< b(x)/2
>> [A,F]=f1('100','0','0', 0, 1, 48, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
ΙΙΙ)κ~=b(x)/2
>> [A,F]=f1('100','0','0', 0, 1, 49, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
ΙV)κ>b(x)/2
>> [A,F]=f1('100','0','0', 0, 1, 50, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
V)κ>> b(x)/2
>> [A,F]=f1('100','0','0', 0, 1, 500, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
B)
Ι)κ<< b(x)/2
>> [A,F]=f1('3628-1987','0','0', 0, 1, 50, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
ΙI)κ< b(x)/2
>> [A,F]=f1('3628-1987','0','0', 0, 1, 800, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
ΙI)κ~= b(x)/2
>> [A,F]=f1('3628-1987','0','0', 0, 1, 821, 1,0,1); >> A\F; >> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) |
ΙV)κ> b(x)/2
>> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) >> [A,F]=f1('3628-1987','0','0', 0, 1, 900, 1,0,1); >> A\F; |
V)κ>> b(x)/2
>> plot(ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf) >> [A,F]=f1('3628-1987','0','0', 0, 1, 10000, 1,0,1); >> A\F; |