ΕΠΙΣΤΗΜΟΝΙΚΟΣ ΥΠΟΛΟΓΙΣΜΟΣ Ι

Ονοματεπώνυμο: Αραβανής Κων/νος  Α .Μ.: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=

·         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+1j   => 

 -1/(h^2) - bj+1/(2*h)= -1/(h^2) + bj/(2*h)  =>

- bj+1 =  bj   =>

 

bj=0

 

Β) Για να είναι θετικά ορισμένο ή ημιορισμένο ένα τέτοιο μητρώο θα πρέπει να ισχύει: x’Ax>0 για κάθε μη μηδενικό διάνυσμα x.

 

Γ) Για να είναι Toeplitz ένα τέτοιο μητρώο θα πρέπει να ισχύει ότι: αj=ct, βj=ct και γj=ct

 

αj= αj+1

βj= βj+1

γj= γj+1

 

2/(h^2) +cj =2/(h^2) +cj+1+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)

 

cj = cj+1

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

γ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+1j  και   γ1n  

Μόνο για το πώτο σκέλος δηλαδή θα πρέπει να ισχύει

 -1/(h^2) - bj+1/(2*h)= -1/(h^2) + bj/(2*h)  =>

- bj+1 =  bj   =>

bj=0

που συνεπάγεται άμεσα και το δεύτερο σκέλος

Επομένος ένα τέτοιο μητρώο είναι συμμετρικό εάν:

bj=0

 

 

 

Β) Για να είναι θετικά ορισμένο ή ημιορισμένο ένα τέτοιο μητρώο θα πρέπει να ισχύει: x’Ax>0 για κάθε μη μηδενικό διάνυσμα x.

 

Γ) Για να είναι Toeplitz ένα τέτοιο μητρώο θα πρέπει να ισχύει ότι: αj=ct, βj=ct και γj=ct

 

αj= αj+1

βj= βj+1

γj= γj+1

 

2/(h^2) +cj =2/(h^2) +cj+1+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)

 

cj = cj+1

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

γj= γj+1

 

2/(h^2) +cj =2/(h^2) +cj+1+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)

 

cj = cj+1

bj = bj+1

 

Ζ) Για να είναι μη αντιστρέψιμο ένα τέτοιο μητρώο πρέπει να υπάρχει τουλάχιστον μία ιδιοτιμή του που να είναι 0.

 

 

Ένας συνοπτικός πίνακας που περιγράφει με ποιές προϋποθέσης το Α μπορεί να είναι συμμετρικό, θετικά ορισμένο ή ημιορισμένο, Toeplitz, αραιό, κυκλοτερές ή μη αντιστρέψιμο φαίνεται παρακάτω:

 

 

 

Dirichlet- Dirichlet

Neumann-Dirichlet

Dirichlet- Neumann

Neumann- Neumann

Περιοδικές

Συμμετρικό

bj=0

bj=0

bj=0

bj=0

bj=0

Θετικά ορισμένο ή ημιορισμένο

x’Ax>0 για κάθε μη μηδενικό διάνυσμα x

x’Ax>0 για κάθε μη μηδενικό διάνυσμα x

x’Ax>0 για κάθε μη μηδενικό διάνυσμα x

x’Ax>0 για κάθε μη μηδενικό διάνυσμα x

x’Ax>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))

 

4) Στο σημείο αυτό εξετάζουμε την συμπεριφορά του δείκτη συνάρτησης κ2(Α) καθώς το n=2^i, i=6:1:10 για μητρώα Α με πολλά από τα παραπάνω χαρακτηριστικά και φυσικά με πολλούς διαφορετικούς συνδιασμούς αυτών. Τα αποτελέσματα παρουσιάζονται σε γραφήματα λογαριθμικής κλίμακας.

 

Για την εκτύπωση των γραφημάτων φτιάξαμε την ακόλουθη συνάρτηση που θεωρεί σαν δεδομένο ότι η συνάρτηση που ψάχνουμε είναι η  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;