Scilab

Scilabin käyttö on hyvin samantyylistä kuin Matlabin, jota monet ovat koulussa ehkä käyttäneet. Matlabin hinnoittelu ei vaan ole kovin käyttäjäystävällistä oppilaitosten ulkopuoliseen käyttöön.

Seuraavassa on pari yksinkertaista esimerkkiä Scilabin käytöstä.

Lineaarinen yhtälöryhmä

Lineaarisen yhtälöryhmän ratkaisu on matriisilaskennan “Hello World”-esimerkki.

1x + 2y + 3z = 4
7x – 3y +2z = -5
1x + 0y – 1z = 2

Kirjoitetaan yhtälöryhmä matriisimuotoon:

Ax = b, joka kerrotaan vasemmalta A:n käänteismatriisilla inv(A):

x = inv(A)*b. Lopuksi tarkastetaan, että tulos täsmää:

Ax – b. Tulokseksi pitäisi saada nollavektori.

Ja sama scilabilla:

--> A = [1 2 3; 7 -3 2; 1 0 -1]
 A  = 

   1.   2.   3.
   7.  -3.   2.
   1.   0.  -1.

--> b = [4; -5; 2]
 b  = 

   4.
  -5.
   2.

--> x = inv(A) * b
 x  = 

   0.9333333
   3.1333333
  -1.0666667

--> A * x - b
 ans  =

   1.776D-15
  -8.882D-16
   0.       

Matka pallon pintaa pitkin

Toinen (hieman monimutkaisempi tai ainakin pitempi) esimerkki vektorilaskennan puolelta.

Kuinka pitkä matka on Hervannan uimahallin pihasta Pariisin Notre Damen aukiolle (teoreettisen) maapallon pintaa pitkin? Ja mihin suuntaan pitäisi lähteä?

Lähtötietoja

Paikat (Google maps)
Piha: 61.453161, 23.842881
Aukio: 48.853581, 2.348070

Ikivanhan metrin määritelmän mukaan matka päiväntasaajalta (pohjois)navalle on kymmenen miljoonaa metriä eli 10 000 km.

Laskenta

Matka

• Määritetään laskentakoordinaatisto, jonka origo on maapallon keskipisteessä, x-akselin suunta on päiväntasaajalle nollapituuspiirin kohdalla ja y-akselin suunta on pohjoisnavalle.
• Halkaistaan maapallo kahtia nollapituuspiirin kohdalta ja määritetään tasolle yksikkövektori, joka osoittaa maapallon keskipisteestä päiväntasaajalle.
• Kierretään em. yksikkövektoria leveyspiirin verran (pohjoiseen). Saadaan suuntavektorin y-koordinaatti.
• Viipaloidaan maapallo sitten em. leveyspiirin kohdalta.
• Kierretään edellisessä kierrossa saatua x-koordinaattia pituuspiirin verran (itään). Saadaan suuntavektorin x- ja z-koordinaatit.
• Sama homma toiselle paikalle ja on saatu kaksi suuntavektoria.
• Suuntavektorien (pituus=1) välisen kulman kosini on vektorien pistetulo. Tästä saadaan niiden välinen kulma.
• Kaaren pituus on säde kertaa kulma (radiaaneina). Ongelman ensimmäinen osa on ratkaistu!

Suunta

• Määritetään tason, jolla ovat pihan ja aukion suuntavektorit, suunta sitä vastaan kohtisuoralla vektorilla.
• Määritetään tason, jolla ovat pihan ja pohjoisnavan suuntavektorit, suunta sitä vastaan kohtisuoralla vektorilla.
• Lasketaan edellisten vektoreiden välinen kulma pistetulon avulla.
• Arvataan meneekö edellä saatu kulma idän vai lännen puolelle. Ongelma on ratkaistu!

Ja sama scilabilla:

// Matka

--> Piha = [61.453161, 23.842881];

--> Piha = Piha * %pi/180 //Kulmat radiaaneiksi
 Piha  = 

   1.07256   0.4161368

--> Aukio = [48.853581, 2.348070];

--> Aukio = Aukio * %pi/180 //Kulmat radiaaneiksi
 Aukio  = 

   0.8526558   0.0409816

--> MaapR = 2 * 10000 / %pi //Maapallon säde
 MaapR  = 

   6366.1977

--> x = [1;0]; //Yksikkövektori

--> Lev = rotate(x,Piha(1))
 Lev  = 

   0.477877 
   0.8784267

--> SPiha = Lev(2); //y-koordinaatti talteen

--> Lev(2) = 0;

--> Pit = rotate(Lev,Piha(2))
 Pit  = 

   0.4370938
   0.1931722

--> SPiha = [Pit(1); SPiha; -Pit(2)] //x- ja z-koordinaatit talteen
 SPiha  = 

   0.4370938
   0.8784267
   -0.1931722

--> Lev = rotate(x,Aukio(1))
 Lev  = 

   0.6579855
   0.7530306

--> SAukio = Lev(2); //y-koordinaatti talteen

--> Lev(2) = 0;

--> Pit = rotate(Lev,Aukio(2))
 Pit  = 

   0.6574331
   0.0269577

--> SAukio = [Pit(1); SAukio; -Pit(2)] //x- ja z-koordinaatit talteen
 SAukio  = 

   0.6574331
   0.7530306
   -0.0269577

--> Pistetulo = sum(SPiha.*SAukio)
 Pistetulo  = 

   0.9540496

--> Kulma = acos(Pistetulo)
 Kulma  = 

   0.3043246

--> Matka = Kulma * MaapR
 Matka  = 

   1937.3904

// Suunta

--> RTulo = cross(SAukio, SPiha) //Ristitulo on kohtisuorassa
 RTulo  = 

  -0.1217842
   0.1152147
   0.2483618

--> Pituus = sqrt(sum(RTulo.*RTulo))
 Pituus  = 

   0.2996488

--> STAukio = RTulo / Pituus //Suunta yksikkövektoriksi
 STAukio  = 

  -0.406423 
   0.3844992
   0.828843 

--> SPohj = [0; 1; 0]; //Pohjoisnavan suunta

--> RTulo = cross(SPohj, SPiha) //Ristitulo on kohtisuorassa
 RTulo  = 

  -0.1931722
   0.       
  -0.4370938

--> Pituus = sqrt(sum(RTulo.*RTulo))
 Pituus  = 

   0.4778770

--> STPohj = RTulo / Pituus //Suunta yksikkövektoriksi
 STPohj  = 

  -0.40423  
   0.       
  -0.9146574

--> Pistetulo = sum(STAukio.*STPohj)
 Pistetulo  = 

  -0.5938190

--> Suunta = acos(Pistetulo)
 Suunta  = 

   2.2065933

--> Suunta = Suunta * 180/%pi //Suunta asteiksi
 Suunta  = 

   126.42849

--> Suunta = 360 - Suunta //Lännen puolella
 Suunta  = 

   233.57151

Taajuusvaste

Lainataan XCos-esimerkissä olevan jousi-massa-vaimennin-järjestelmän siirtofunktiota, mutta hieman pienemmällä vaimennuksella, ja määritellään sille taajuusvaste (Bode).

Siirtofunktio: 1 / (s² + s + 20)

Scilab:

--> s = %s
 s  = 

  s

--> G = 1/(s*s+s+20)
 G  = 

              
       1      
   ---------  
   20 +s +s²  

--> sys = syslin('c',G)
 sys  = 

              
       1      
   ---------  
   20 +s +s²  

--> bode(sys,0.1,5)

Ja tulos:

x