Jakobiyen matrisi ve 3. derece polinom

[Tesla Coil]

Soruna binaen aşağıdaki dökümanı hazırladım umarım faydalı olur. Adımları tekrarlıyarak istediğin kadar iterasyon yapabilirsin ben sadece 1. iterasyon için yaptım. Başarılar

https://www.dropbox.com/s/y3gdvl7z720bfya/3.mertebejacop.pdf

http://s3.dosya.tc/server23/dpzAMk/3.mertebejacop.pdf.html




Kaynaklar
[1] F. M. Nuroğlu Power System Analysis ders notları,
[2] Power System Analysis – P.S.R. Murty

[Cemre.]

Hocam uğraşmışsınız gerçekten çok teşekkürler...

Ben de şöyle birşeye ulaştım sonuç olarak, ancak sonuçtan bayağı bir uzaklaşıyorum.. Örneği bir kontrol edebilecek olan var mı?

program saycoz;
uses crt;
type MATRIS = array [1..20,1..20] of real;
Var A0,A1,A2,A3,X1b,X2b,X3b,f1,f2,f3:real ;
n,i,j,k,q,u,M:integer;
W,Xy,Xc,B:MATRIS;
p,t,det:real;
Label 1;



begin
Writeln('Hesaplanacak polinom: A0x^3 + A1x^2 + A2x + A3');
write('A0=');
readLn(A0);
write('A1=');
readLn(A1);
write('A2=');
readLn(A2);
write('A3=');
readLn(A3);

A1:=A1/A0;
A2:=A2/A0;
A3:=A3/A0;

writeln('Baslangic degerleri,');
write('X1=');
readLn(X1b);
write('X2=');
readLn(X2b);
write('X3=');
readLn(X3b);
writeln('Kaç iterasyon yapılacak?');
readLn(u);

for q:=1 to u do
begin

f1:=A1+X1b+X2b+X3b;
f2:=A2-(X1b*X2b)-(X1b*X3b)-(X2b*X3b);
f3:=A3+(X1b*X2b*X3b);

W[1,1]:=1; 
W[1,2]:=1;
W[1,3]:=1;
W[2,1]:=-X2b-X3b;
W[2,2]:=-X1b-X3b;
W[2,3]:=-X2b-X1b;
W[3,1]:=X2b*X3b;
W[3,2]:=X1b*X3b;
W[3,3]:=X2b*X1b;

N:=3;M:=3;det:=1;
B[1,1]:=1; B[1,2]:=0; B[1,3]:=0;
B[2,1]:=0; B[2,2]:=1; B[2,3]:=0;
B[3,1]:=0; B[3,2]:=0; B[3,3]:=1;
 

 for q:=1 to 3 do begin  //ana döngü
 
 for j:=1 to 3 do begin  //inverse baslangıc
 B[j,j]:=1; 
 t:=W[j,j];
 det:=det*t; 
             for k:=1 to N do begin
             W[j,k]:=W[j,k]/t;
             B[j,k]:=B[j,k]/t;
             end;
             for i:= 1 to M do begin
             p:=W[i,j];
             if i=j then goto 1;
             for k:=1 to 3 do begin
             B[i,k]:=B[i,k]-p*B[j,k];
             W[i,k]:=W[i,k]-p*W[j,k];
             end;
1:end;
end;
for i:=1 to 3 do begin
for j:=1 to 3 do writeln('inv',i,j,' ',B[i,j]:8:4);
end;
//inverse sonu

Xc[1,1]:=(B[1,1]*f1)+(B[1,2]*f2)+(B[1,3]*f3);
Xc[2,1]:=(B[2,1]*f1)+(B[2,2]*f2)+(B[2,3]*f3);
Xc[3,1]:=(B[3,1]*f1)+(B[3,2]*f2)+(B[3,3]*f3);

Xy[1,1]:=((X1b)-(Xc[1,1]));
Xy[2,1]:=((X2b)-(Xc[2,1]));
Xy[3,1]:=((X3b)-(Xc[3,1]));

X1b:=Xy[1,1];
X2b:=Xy[2,1];
X3b:=Xy[3,1];
writeLn(q,'inci iterasyon sonucu: x1=',X1b);
writeLn(q,'inci iterasyon sonucu: x2=',X2b);
writeLn(q,'inci iterasyon sonucu: x3=',X3b);
end;
end;
readln;
end.

polinomu A0=1 A1=-6 A2=11 A3=-6 giriyoruz, başlangıç değerleri X1=0.5 X2=1.5 X3=2.5 olacak şekilde, zaten köklerin de x1=1,x2=2,x3=3 olarak gelmesi gerekiyor ancak bizim değerler 3üncü iterasyona kadar hesapladığımda 20küsürlere fırlıyor.. Gerçekten beynim dondu hata bulamıyorum kim bilir kaç kez kontrol ettim

Ha bir de bu 3iterasyondan fazlasında sapıtıyor, bu pascal ile ilgili bir durum olabilir mi acaba?

[Tagli]

Daha önce sözü geçen bulanık noktaları netleştirdin mi?

1) Kaç adet fonksiyon olacak, 1 mi 3 mü? 1x3 matrisin tersini normal yöntemle alamazsın.
2) J analitik türev alınarak mı hesaplanacak yoksa sayısal türev alınarak mı?

Çözümün, başlangıç değerlerine bağlı olarak yakınsamama ihtimali de var. Eğer birden fazla çözüm varsa, hangisine yakınsayacağı yine başlangıç değerlerine bağlı.

Ayrıca, genelde iterasyon sayısı önceden verilmez. Sadece bir üst sınır belirlenir. İterasyonlarda bulunan değerler birbirine yakın çıkmaya başladıklarında işlem kesilir.

[Tesla Coil]

Yolladığım dosyayı inceledin mi? Pascalla daha önce hiç uğraşmadım ama değişken tanımlamaların çok karmaşık geldi.

Program tam olarak ne yapıyor anlamadım. Matematik konusunda anlamadığın bir nokta varsa yardımcı olayım. Eğer epsilon gibi bir hata opsiyonu tanımlıyacaksan newton raphsonda genelde ∆x0 olarak tanımlanan hata değerini kullanabilirsin ki genelde ona bakarlar. ∆x0 ve ∆x1 arasındaki farka bakılır diye biliyorum 2 iterasyon arasındaki farka değil.

[Cemre.]

Hocam fonksiyon sayısı 3, ben yanlış anlamışım zaten direk polinomun değil, programda da hesaplandığı gibi f1 f2 f3'ün kısmi türevleri alınacak, bu noktada sıkıntı yok. Türevleri analitik yolla hesapladık zaten programda da görülüyor. Başlangıç değerlerini hiç değiştirmedim o konuda haklı olabilirsiniz hemen deniyorum.. Aslında bir epsilon değeri belirleyip iki iterasyon sonucu arasındaki farkın epsilondan büyük küçük kontrolü yaparak döngüden çıkartmamız gerekir ama, böyle bir şey istenmiyor bildiğim kadarıyla.
mesaj birleştirme:: 26 Mayıs 2014, 23:40:34
Evet hocam değişkenler biraz karmaşık, arkadaşın tanımladığı yerden devam ettirmemden kaynaklı ben de baya bir zorluk çektim. Ancak matematiksel işlemler düşünülürse neyin ne olduğu anlaşılıyor. Ben direk sorayım daha kolay olur sanırım.

for j:=1 to 3 do begin  //inverse baslangıc
 B[j,j]:=1; 
 t:=W[j,j];
 det:=det*t; 
             for k:=1 to N do begin
             W[j,k]:=W[j,k]/t;
             B[j,k]:=B[j,k]/t;
             end;
             for i:= 1 to M do begin
             p:=W[i,j];
             if i=j then goto 1;
             for k:=1 to 3 do begin
             B[i,k]:=B[i,k]-p*B[j,k];
             W[i,k]:=W[i,k]-p*W[j,k];
             end;
1:end;
end;
for i:=1 to 3 do begin
for j:=1 to 3 do writeln('inv',i,j,' ',B[i,j]:8:4);
end;
//inverse sonu

matris tersi alan kod parçacığı, acaba hatalı mı değil mi? (B matrisi birim matris, W matrisi jakobiyen matris)

Xc[1,1]:=(B[1,1]*f1)+(B[1,2]*f2)+(B[1,3]*f3);
Xc[2,1]:=(B[2,1]*f1)+(B[2,2]*f2)+(B[2,3]*f3);
Xc[3,1]:=(B[3,1]*f1)+(B[3,2]*f2)+(B[3,3]*f3);

Xy[1,1]:=((X1b)-(Xc[1,1]));
Xy[2,1]:=((X2b)-(Xc[2,1]));
Xy[3,1]:=((X3b)-(Xc[3,1]));

burada, B yani Jakobiyen matrisin tersi, F matrisiyle çarpılıyor daha sonra başlangıç değerinden tek tek çıkarılıyor ve çıkan sonuç, Xy sonuç değişkenine atanıyor, o da zaten her iterasyonun sonucu ve bu sonuc Xb değerlerine yani başlangıç değerlerine tekrar atanıp başa dönülüyor böylece bir iterasyon tamamlanmış oluyor.. Benim gördüğüm matematik kısmın gösterdiğiniz kısımla aynı olduğu ama bir problem var ki sonuç alakasız çıkıyor :/

--------
Hocam bir de başlangıç değerlerini 1,2,3 olarak verirsem yani kökleri direk yazarsam her iterasyonun sonucu 1,2,3 olarak geliyor, yani sanki matematik kısmında da bir sorun yok, acaba sayısal türev aldırarak mı çözmeli jakobiyen matrisi?

[Tesla Coil]

Ayrıca 1. iterasyonda bulduğun değer 2. iterasyonunun başlangıç değeri.

Her yeni iterasyonda en son bulduğun sonucu başlangıç değer olarak kullanmalısın. İlk girilen başlangıç değerlerin sadece 1. iterasyon için 2. iterasyondaki başlangıç değerlerin 1. iterasyonun sonucları olması gerekiyor buna dikkat ediyomusun?
x1= x0+∆x0
x2=x1+∆x1 diye devam ediyor.


mesajı 2 saat önce girdim hala moderatör onaylicak.!!

edit2: Mesajı kaç saat önce girdiğimi hatırlamıyorum. Moderatör beyfendi onaylamadı. Özelden yazıyım dedim sanırsam böyle bir yetkimde yokmuş.

Gelelim konumuza;

J matrisin her iterasyonda değişecek.
1. iterasyonda başlangıç değerlerin var ve bunu j de kullanıyosun
2. iterasyonda 1. iterasyonun sonuçlarını kullanırsın.

[Tesla Coil]

Ayrıca 1. iterasyonda bulduğun değer 2. iterasyonunun başlangıç değeri.

Her yeni iterasyonda en son bulduğun sonucu başlangıç değer olarak kullanmalısın. Başlangıç değerlerin sadece 1. iterasyon için 2. iterasyondaki başlangıç değerlerin 1. iterasyonun sonucları olması gerekiyor buna dikkat ediyomusun??
x1= x0+∆x0
x2=x1+∆x1 diye devam ediyor.

[Tagli]

Pascal bilmediğimden olacak herhalde, kodu çözemedim. Zaten biraz da karışık yazılmış. Ben C'de matris tersi ve determinant almak için epey kod yazdığımı hatırlıyorum. Ama kodun uzun olması, büyük matrisler için de çalışabilecek genel bir kod yazmış olmamdan kaynaklanıyordu. 3x3 için kestirme yollar olabilir. Belki sendeki kod onlardan birini kullanıyordur. Ama bence, J'yi ve tersini ekrana bastır ve başka bir yazılımla tersinin doğru alınıp alınmadığını kontrol et. Daha sonra da matris çarpma ve toplama işlemlerini kontrol et.

[Cemre.]

Dediğiniz gibi kontrol ettim hocam inverste sorun yokmuş..
Benim bir sorum olacak bu arada, şu F matrisi, yani Tesla hocamın dropbox'a attığı dökümandaki deltaC matrisi programın en başında verdiğimiz başlangıç değerleriyle hesaplanıp hep sabit mi kalacak yoksa her iterasyondan sonra aldığımız başlangıç değerleriyle tekrar mı hesaplanacak?
mesaj birleştirme:: 27 Mayıs 2014, 00:36:04
Sonunda oldu sanırım.. Hocam çok teşekkürler emekleriniz için, Tesla hocam hazırladığınız döküman için ayrıca teşekkürler..

Bu J matrisi beni öldürecekti bu gece.. Meğer J'yi bir kez hesaplayıp birdaha dokunmamak gerekiyormuş yada ben bayaa büyük bir hatanın içerisindeyim ancak ilk kez tam olarak kökleri buldurabildim ve başlangıç değerlerimi değiştirdiğimde de hiç şaşmadan buluyor köklerimi, bilmiyorum umarım son hali olur bu, gerçekten yordu çünkü.. Son kez bakar mısınız hocam doğru mudur J sabit mi kalacak??

program saycoz;
uses crt;
type MATRIS = array [1..20,1..20] of real;
Var A0,A1,A2,A3,X1b,X2b,X3b,f1,f2,f3:real ;
i,j,k,q,u:integer;
W,Xy,Xc,B:MATRIS;
p,t,det:real;
Label 1;



begin
Writeln('Hesaplanacak polinom: A0x^3 + A1x^2 + A2x + A3');
write('A0=');
readLn(A0);
write('A1=');
readLn(A1);
write('A2=');
readLn(A2);
write('A3=');
readLn(A3);

A1:=A1/A0;
A2:=A2/A0;
A3:=A3/A0;

writeln('Baslangic degerleri,');
write('X1=');
readLn(X1b);
write('X2=');
readLn(X2b);
write('X3=');
readLn(X3b);
writeln('Kaç iterasyon yapılacak?');
readLn(u);

det:=1;
B[1,1]:=1; B[1,2]:=0; B[1,3]:=0;
B[2,1]:=0; B[2,2]:=1; B[2,3]:=0;
B[3,1]:=0; B[3,2]:=0; B[3,3]:=1;

W[1,1]:=1; 
W[1,2]:=1;
W[1,3]:=1;
W[2,1]:=-X2b-X3b;
W[2,2]:=-X1b-X3b;
W[2,3]:=-X2b-X1b;
W[3,1]:=X2b*X3b;
W[3,2]:=X1b*X3b;
W[3,3]:=X2b*X1b;

 for q:=1 to u do begin  //ana döngü

f1:=A1+X1b+X2b+X3b;
f2:=A2-(X1b*X2b)-(X1b*X3b)-(X2b*X3b);
f3:=A3+(X1b*X2b*X3b);

 for j:=1 to 3 do begin  //inverse baslangıc
 t:=W[j,j];
 det:=det*t; 
             for k:=1 to 3 do begin
             W[j,k]:=W[j,k]/t;
             B[j,k]:=B[j,k]/t;
             end;
             for i:= 1 to 3 do begin
             p:=W[i,j];
             if i=j then goto 1;
             for k:=1 to 3 do begin
             B[i,k]:=B[i,k]-p*B[j,k];
             W[i,k]:=W[i,k]-p*W[j,k];
             end;
1:end;
end;
for i:=1 to 3 do begin
for j:=1 to 3 do writeln('inv',i,j,' ',B[i,j]:8:4);
end;
//inverse sonu

Xc[1,1]:=(B[1,1]*f1)+(B[1,2]*f2)+(B[1,3]*f3);
Xc[2,1]:=(B[2,1]*f1)+(B[2,2]*f2)+(B[2,3]*f3);
Xc[3,1]:=(B[3,1]*f1)+(B[3,2]*f2)+(B[3,3]*f3);

Xy[1,1]:=((X1b)-(Xc[1,1]));
Xy[2,1]:=((X2b)-(Xc[2,1]));
Xy[3,1]:=((X3b)-(Xc[3,1]));

X1b:=Xy[1,1];
X2b:=Xy[2,1];
X3b:=Xy[3,1];
writeLn(q,'inci iterasyon sonucu: x1=',X1b:5:3);
writeLn(q,'inci iterasyon sonucu: x2=',X2b:5:3);
writeLn(q,'inci iterasyon sonucu: x3=',X3b:5:3);
end;
readln;
end.

[Tagli]

Hem J'nin hem de F'nin her adımda, o adımın tahmini değişkenleri (x'ler) için yeniden hesaplanması gerekir.

[Cemre.]

Evet hocam formüle bakıldığında öyle, ancak dediğiniz gibi yaptığımda olmuyor da olmuyor, böyle nedense oldu
mesaj birleştirme:: 27 Mayıs 2014, 00:58:22
Hocam şuan hatayı farkettim, dediğiniz gibi oluyor, ben matrisin tersini ikinci kez hesaplarken birim matrisi yeniden B'ye eşitlemiyormuşum, şimdi oldu işte.. Tekrardan çok teşekkürler
mesaj birleştirme:: 27 Mayıs 2014, 01:09:46
Bu arada, son hali, belki birinin işine yarar...

program saycoz;
uses crt;
type MATRIS = array [1..20,1..20] of real;
Var A0,A1,A2,A3,X1b,X2b,X3b,f1,f2,f3:real ;
i,j,k,q,u:integer;
W,Xy,Xc,B:MATRIS;
p,t,det:real;
Label 1;



begin
Writeln('Hesaplanacak polinom: A0x^3 + A1x^2 + A2x + A3');
write('A0=');
readLn(A0);
write('A1=');
readLn(A1);
write('A2=');
readLn(A2);
write('A3=');
readLn(A3);

A1:=A1/A0;
A2:=A2/A0;
A3:=A3/A0;

writeln('Baslangic degerleri,');
write('X1=');
readLn(X1b);
write('X2=');
readLn(X2b);
write('X3=');
readLn(X3b);
writeln('Kaç iterasyon yapılacak?');
readLn(u);


for q:=1 to u do begin  //ana döngü (iterasyon işlemi)

det:=1;
B[1,1]:=1; B[1,2]:=0; B[1,3]:=0;
B[2,1]:=0; B[2,2]:=1; B[2,3]:=0;
B[3,1]:=0; B[3,2]:=0; B[3,3]:=1;

W[1,1]:=1; 
W[1,2]:=1;
W[1,3]:=1;
W[2,1]:=-X2b-X3b;
W[2,2]:=-X1b-X3b;
W[2,3]:=-X2b-X1b;
W[3,1]:=X2b*X3b;
W[3,2]:=X1b*X3b;
W[3,3]:=X2b*X1b;


for j:=1 to 3 do begin  //inverse baslangıc
t:=W[j,j];
det:=det*t;
   for k:=1 to 3 do begin
   W[j,k]:=W[j,k]/t;
   B[j,k]:=B[j,k]/t;
   end;
   for i:= 1 to 3 do begin
   p:=W[i,j];
   if i=j then goto 1;
   for k:=1 to 3 do begin
   B[i,k]:=B[i,k]-p*B[j,k];
   W[i,k]:=W[i,k]-p*W[j,k];
   end;
1:end;
end;
writeLn('Jakobiyen matrisin tersi;');
for i:=1 to 3 do begin
for j:=1 to 3 do writeln(i,',',j,' ',B[i,j]:8:4);
end;
writeLn();
//inverse sonu

f1:=A1+X1b+X2b+X3b;
f2:=A2-(X1b*X2b)-(X1b*X3b)-(X2b*X3b);
f3:=A3+(X1b*X2b*X3b);

Xc[1,1]:=(B[1,1]*f1)+(B[1,2]*f2)+(B[1,3]*f3);
Xc[2,1]:=(B[2,1]*f1)+(B[2,2]*f2)+(B[2,3]*f3);
Xc[3,1]:=(B[3,1]*f1)+(B[3,2]*f2)+(B[3,3]*f3);

Xy[1,1]:=((X1b)-(Xc[1,1]));
Xy[2,1]:=((X2b)-(Xc[2,1]));
Xy[3,1]:=((X3b)-(Xc[3,1]));

X1b:=Xy[1,1];
X2b:=Xy[2,1];
X3b:=Xy[3,1];

writeLn(q,'.iterasyon','  X1=',X1b:5:3,'  X2=',X2b:5:3,'  X3=',X3b:5:3);
writeLn();
end;
readln;
end.

[Tesla Coil]

Dostum "acemi üye" sıfatında olduğum için yazdığım cevaplar gelmedi. Sanırım sen sorunlarını çözmüşsün pürüz var mı ?

[Cemre.]

Şuan bir sıkıntı yok hocam, çok sağ olun tekrardan, sınıfın yarısı nasiplendi