Matlab-harjoitus, harj0

Harj0: Matriisilaskennan ja Matlabin kertaus. Heti 1. tiistaina alkaa. Lue rinnakkain vastaavia Maple-operaatioita [HAM]-kirjan Luvusta 7: Matriisit ja vektorit s. 143 ->

Diff. yhtälöt ja lineaarialgebra, LAODE

Kts ./matlab/laode.m Itseopiskelua ja kertausta

Luku 1

1.2

sivu 6

format compact
A=[2 3 1;1 4 7]
A(1,:),A(2,:)
A(:,3)
A(1:2,2:3)

A([1 1 2 1 2 2],1:2)

Tehtävä 1, s. 6 Muodosta 3x4-matriisi A ja poimi joitakin alkioita, rivien lineaarikombinaatioita ...

Loput tehtävät tässä kohdassa liian lapsellisia.

1.3 Special kinds of matrices

diag(1:3)
eye(5)
zeros(4,5)
ans'

1.4 The geometry of vector operations

x=[1 2];y=[-2 3];addvec(x,y)
x=[1, 0, 2]; y=[-1, 4, 1]; % Pilkkujen käyttö on yhdenmukaisempaa Maplen kanssa
addvec3(x,y); grid

addvec on laode-funktio. Vertailun vuoksi kannattaa katsoa [HAM] s. 145.

Vektorin piirto pelkkänä janana on tietysti äärimmäisen yksinkertainen asia:

x=1:2;
plot([0 x(1)],[0,x(2)])
x=1:3;
plot3([0,x(1)],[0,x(2)],[0,x(3)]);grid

Monikäyttöinen vecplot on peräisin kirjasta [RiSm]. Se on matlab-hakemistoissamme

help vecplot

  a -- 2 x n matriisi: Sarakkeet vektorien loppupisteet.
  v -- 2 x 1 sarakevekt: yhteinen alkupiste (opetus: Origo).
 
  Esim: 
  t=pi/3;R=[cos(t),-sin(t);sin(t),cos(t)];v=[.5;1];vecplot([v,R*v])
  axis([-1 1 0 1]);
  
  clf;v=[1;0]
  w=v;w=R*w;vecplot([v,w]);axis([-1 1 -1 1]);hold on
  w=R*w;vecplot([v,w]); Iteroidaan tätä riviä.

Normi Matlabissa norm:n oletuksena on 2-normi, eli euklidinen normi

for n=1:10
   norm(ones(n,1))
end

Maplessa oletusnormina on "ääretön-normi". Euklidinen saadaan näin:

norm(vector([a,b,c]),2); Kts. [HAM] s. 146 Matlabissa toimivat:


>> norm((1:10),inf)     
ans =
    10
>> norm((1:10),2)  
ans =
   19.6214
>> norm((1:10))  
ans =
   19.6214

Sisätulo (pistetulo) dot, ristitulo cross

x=[1 4 2];y=[2 3 -1];
dot(x,y)
x*y'       % myös näin, mieti miksi!
sum(x.*y)  % ja näinkin, mieti taas! 
cross(x,y)

Tehtävä Määritä vektorin u kohtisuora projektio vektorilla v. Testaa nyt vaikka näillä:

u=1:3
v=[1,0,0];% v=[0,1,0] % v=[0,0,1]

Luku 2, Lin. equ s. 23

Ratkaistaan joitakin systeemejä Matlabilla Tehtäviä s. 23

8. Seuraavat komennot lataavat matriisin A ja pystyvektorin b. Kyseessä ovat laode-polun varrella olevat komentotiedostot e2_1_4.m ja e2_1_5.m . Katso niiden sisältö suorittamalla esim.
type e2_1_4
```
e2_1_4 
e2_1_5 
```
Ratkaise lineaarinen yhtälösysteemi ja tarkista kertomalla. Laske erityisesti erotus A*x - b . format long . Laske myös erotuksen normi.
9..12 Matriisin alkioita muutellaan. Voit kokeilla jotain lyhyesti (älä jää pitkäksi aikaa).

Muutetaan matriisin alkio (2,3) niin, että joudutaan pulmatilanteeseen.

A(2,3)=1
x=A\b

Laske A:n determinantti ( det ) ja "häiriöalttius" ( cond ).

Suoritustehtäviä, aihe: lin. yhtsyst. soveltamista

Tässä muutama pikku sovellus, itse ratkaiseminen otetaan mustana laatikkona:

             x=A\b

Määritä 3. asteen polynomi, joka toteuttaa seuraavat interpolaatioehdot:
```
     p(1)=2, p(2)=3, p'(-1)=-1, p'(3)=1.
```
Piirrä polynomi sopivalla välillä ja samaan kuvaan myös "datarinkulat".
Kyseessä on ns. Hermiten interpolaatio , jossa asetetaan ehtoja niin funktion arvoille kuin derivaatoille. (Välimuoto Taylorin ja Lagrangen suhteen.)
Ohje: Kirjoita kynällä yhtälösysteemi, syötä sen matriisi Matlabiin ja ratkaise. Polynomin arvojahan lasketaan polyval-funktiolla. Täältä vaikka lisää tai help polyval .

Muodosta lineaarinen yhtälösyteemi alla olevalle virtapiirille ja ratkaise virrat.

                         
               8 V
             ----->              Kirchhofin virtalaki:
        I1    |              
    ----<---- | |---------   Virtapiirin jokaisessa pisteessä pätee:
    |         |          |   tulevien virtojen summa =  lähtevien 
    <                    <   summa.
    > 1 Ohmi      1 Ohmi >
    <                    <       Kirchhofin jännitelaki:
    >                    >
    |      0.5 Ohmia     |   Virtapiirin jokaisessa suljetussa silmukassa
    ------\/\/\/--->-----|   ulkoisen jännitelähteen jännite =
    |              I3    |   jännitehäviöiden summa.
    <                    |
    >                    |      Ohmin laki: U=RI
    < 2 Ohmia         I2 \/ 
    >                    |
    |         |          |
    ----<---- | |---------
              |          
            16 V
           ------>

Tasoalueen reunoilla on oheisen kuvan mukaiset lämpötilat.
```
          50   20   10   0
     |----|----|----|----|----|
     |                        |
  80 |   T1   T2   T3   T4    | 80
     |                        |
  60 |   T5   T6   T7   T8    | 60
     |                        |
  40 |   T9   T10  T11  T12   | 40
     |                        |
  20 |   T13  T14  T15  T16   | 20
     |                        |
     |----|----|----|----|----|
          0    0    0    0
```
Oletetaan, että sisäpisteiden lämpötilat T1,...,T16 saadaan naapuripisteiden lämpötilojen keskiarvona. (Kullakin sisäpisteellä on 4 naapuria, pohjois-, etelä-, itä-, länsi.) Muodosta lineaarinen yhtälösysteemi lämpötilojen T1,...,T16 ratkaisemiseksi. Rakenna yhtälösysteemi itsellesi ensin kynää ja paperia käyttäen ja syötä se sitten Matlabiin. Käytä ratkaisuun Matlabin valmista ratkaisijaa T=A\b . Piirrä lämpötilafunktion kuvaaja. Tämä käy muotoilemalla ratkaisuvektori T 4x4 - matriisiksi TMAT . ja soveltamalla surf -tyyppistä funktiota. Huom: TMAT -matriisin saa kätevästi komennon reshape avulla. Voit vielä liikutella pintaasi vaikkapa tähän tapaan:
```
for j=1:10;view(-20,10*j),pause,end;  
```
Tosin pyörittely käy nykyisin hyvin myös valikosta hiiren avulla. Tehtävä on tarkoitus tehdä Matlabilla. Maple-ohje vain asian lisäharrastukseen.
Huom! Tässä on kyse ihan oikeasta numerisesta menetelmästä ns. Laplacen osittaisdifferentiaaliyhtälön ratkaisemiksi, ns. differenssimenetelmästä.
Maple-ohje Systeemin voi ratkaista yhtä hyvin käyttäen Maple:n linalg[linsolve] - komentoa, jolloin muotoilu (vrt reshape edellä) voidaan tehdä ihan vain arkisella matrix komennolla. Piirtoon plots[matrixplot]

2.2 The geom. of low-dim. sols s. 24 ..

Suorien leikkauspiste:

x+y=7
-x+3y=1 Tässä komennot, selvitä, leikkaa/liimaa Matlabiin.

x=[-3,7];  % Ei tarvita 100:aa pistettä, kaksi riittää.
y=7-x
plot(x,y,'b')
xlabel('x');ylabel('y');
hold on
y=(1+x)/3;
plot(x,y,'r')
axis('equal')
grid

Jatketaan vielä laskemalla suorien leikkauspiste ja piirtämällä se rinkulalla.

A=[1 1;-1 3];b=[7;1];
x=A\b
plot(x(1),x(2),'o')
hold off

Tästä voitaisiin jatkaa tasoihin, saataisiin samalla meshgrid-kertaus (ss. 28-30). Ehkä jätetään se väliin ja katsotaan meshgrid-asiaa tuon tuostakin muuten. Tässä (s. 30) kerrataan myös yleinen piirto "nonlinear functions". Lue kirjasta ja kokeile jos tunnet tarvetta. Linkkejä ... Tehtäviä ss 32-33

7. Ratkaise yhtälöpari
```
    x + 4y = -4
   4x + 3y = 4
```
sekä graafiseti että numeerisesti. Kokeile zoomia. Piirrä numeerisen ratkaisun perusteella rinkula oikeaan paikkaan.
11. Piirrä funktion f(x)=2-x*sin(x^2-1) kuvaaja välillä [-2,3]. Kuinka monta maksimia tällä välillä? Lapsellinen kysymys, otetaan vähän "miehekkäämpi": Määritä maksimi (100 pisteen resoluutiolla) ja laske vastaava x:n arvo. No, näitä oli V2-kevään alussa. y-vektorin suurin löytyy max-funktiolla, olkoon se muuttujassa M. Näppärintä siihen max-kohtaan lienee tämä: xmax=x(y==M). (Indeksoidaan bittivektorilla, jolloin poimitaan se (ne) joissa bitti on päällä .)
Huom! Muistathan ne pisteet: .* ja .^ Lisäksi kärsimystä voi minimoida puolipisteellä, kun joudut lakemaan 100 (tai 200) numeroa, joita et pala halusta nähdä. Piirrä se maksimipiste taas rinkulalla samaan kuvaan.

2.3 Gauss

Rivioperaatiot:


A(i,:)=c*A(i,:);
A(j,:)=A(j,:)+c*A(i,:);
A([i j])=A([j i]);

e2_3_8

Matlab-funktioita

Tehtävä (ehkä vähän myöhemmin) Värähtelyliikkeen yhteydessä tarvitsemme usein kaavaa

      a sin wt + b cos wt = C cos(wt -delta)

Tee Matlab-funktio, joka laskee amplitudin C ja vaihekulman delta, kun syötteenä annetaan kertoimet a ja b. Jos olet unohtanut cos:n yhteenlaskukaavan, voit kysyä maplelta

cos(omega*t-delta);%=expand(%);

Huomaa, että vaihekulma ei tule yksikäsitteisesti määrätyksi lausekkeesta arctan(b/a), vaan koordinaattineljännes pitää päätellä a:n ja b:n merkeistä. Tätä varten on Matlabissa funktio atan2. (Myös angle(a+i*b) hoitaa homman).

%%%%%%%%%%%
function[C,delta,lauseke] = amplitmuoto(a,b,w,t)
% Kutsu 1: [C,delta,lauseke] = amplitmuoto(a,b,w,t)
% Kutsu 2: [C,delta] = amplitmuoto(a,b)
% Syöte: a --  cos-termin kerroin
%        b --  sin-termin kerroin
%        w --  kulmataajuus (valinnainen)
%        t --  vektori, jossa lausekkeen arvot lasketaan (valinnainen)
C= ... ;
delta = ... ;
if nargout < 3
   return;

if nargin > 2
   lauseke = ...
end;

3.7

Tehtäviä s. 101

9
```
A:=matrix(3,3,[1,a,b,0,1,c,0,0,1]);
AI:=augment(A,diag(1,1,1));
gaussjord(AI);
```
a) Onko matriisi kääntyvä millä tahansa arvoilla a,b,c ?
b) Määritä käänteismatriisi niillä joilla on (kenties kaikilla?) Vastaa a)- kysymykseen ilman minkään ohjelman apua. b)-kohtaan kannattaa käyttää Maplea.

A=[2 1 3;1 2 3;5 1 0]
AI=[A eye(3)]
IA=rref(AI)
B=IA(:,4:6)
[A*B,B*A]  % tarkistus

C=[1 0 3;-1 2 -2;0 2 1]
inv(C)
CI=[C,eye(3)]
format rational; format compact
rref(CI)
» rref(CI)
ans =
      1            0            3            0           -1            1      
      0            1           1/2           0            0           1/2     
      0            0            0            1            1           -1

Rref-muodon alkio (3,3) = 0. Vastaavat bmato-alkiot eivät ole 0, joten ratkaisuja ei ole.

> A:=matrix(3,3,[1,epsilon,3,-1,2,-2,0,2,1]);
> AI:=augment(A,diag(1,1,1));
> gaussjord(AI);

Teht. 3.8. s. 105

11--14

A=[0 -2;2 0];det(A)
4
A=[-0.5 -0.5;0.7 0.7];det(A)
0
A=[-1 -0.5;-2 -1],det(A)  % toinen näistä riittää.
0
format short
A=[0.7071 0.7071;-0.7071 0.7071],det(A)

Lineaarikuvauksia

Liittyy luentoaiheeseen to 14.9.00


%function [A,Theta,c]=kiertod(a,b)
%A=[a -b;b a];
%Theta=atan2(b,a);
%c=sqrt(a^2+b^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

format compact
[A,T,c]=kiertod(3,4)
A =
     3    -4
     4     3
T =
    0.9273
c =
     5
     
T/pi*180
ans =
   53.1301     
     
   map
   Tarkistetaan, lähdetään hyvin lyhyestä liikkeelle (0.1,0).

Lineaarikuvauksen havainnollistus: yksikköympyrän kuva

Rakennetaan maanläheinen skripti yksikköympyrän kuvautumisesta.

   clf
   A=diag(1:2)              % Vaihda A siten kuin haluat.
   t=linspace(0,2*pi,50);   % Kulmat
   ympyra=[cos(t);sin(t)];  % 1-ymp. pisteet: 1. rivi: x-pisteet
                            %                 2. rivi: y-pisteet
   kuva=A*ympyra;           % Kuvapisteet A:lla kerrottaessa

   plot(ympyra(1,:),ympyra(2,:),'r')  % Ympyrä punaisella 'r'
   hold on  

   plot(kuva(1,:),kuva(2,:),'b');     % Kuva sinisellä
   axis('square');axis([-4 4 -4 4])
   grid;shg;
   disp('Valitse punaiselta ympyrältä piste, jonka kuvan haluat nähdä')
   [x,y]=ginput(1)                    
plot([0,x],[0,y],'r')
uv=A*[x;y];
plot([0,uv(1)],[0,uv(2)],'b')

Saadaan tällainen:

[linkuv1.gif]

Tähän olis kiva laittaa for-silmukka, joka antaisi mahdollisuuden valita useita pisteitä ja katsoa niiden kuvautumista