Laboration i Matlab. Uppgift 1. Beskrivning

(1)

Laboration i Matlab

Uppgifterna i denna laboration kan innehålla fysik och matematik som ni inte känner till, men det kommer ni inte att behöva för att kunna lösa uppgifterna. Uppgifterna är skrivna så att ni får en del hjälp på vägen, men vissa delar får ni göra själva. Det går mycket lättare om in har läst kompendiet Introduktion till Matlab innan ni genomför laborationen.

Uppgift 1

Beskrivning

Matlab är mycket effektiv på att lösa stora ekvationssystem med flera obekanta, och det finns många fysikaliska problem där man behöver detta. Till exempel inom elläran så får man ofta stora ekvationssystem av obekanta strömmar och spänningar, som måste lösas. I den första uppgiften ska ni därför lösa ut strömmar och spänningar i ett sådant ekvationssystem, men ni kommer att få hjälp med att ställa upp systemet.

Betrakta det elektriska schemat nedan. Alla strömmar och spänningar finns redan angivna, men det finns flera steg att gå igenom innan ekvationssystemet kan lösas.

Först måste vi använda lite kunskaper i ellära.

Fig. 1. Ett elektriskt schemat med obekanta strömmar och spänningar.

Schemat i Fig. 1 går att räkna ut för hand utan, men nu ska vi illustrera hur det går att beräkna denna krets med Matlab. Till och börja med måste vi identifiera obekanta och inparametrar. I detta fall är spänningarna, V1, V2, V3 och V4, insignaler till kretsen, så dessa måste vi ange ett värde på för att kunna beräkna de övriga (I1-I7, I+, I- och Vut, V+, V-), som då är obekanta. Vidare så måste vi också känna till alla värden på resistanserna, R1-R7, samt förstärkningen A (se nedan).

Vi börjar med en förenkling: Den trekantiga symbolen illustrerar en operationsförstärkare, som förstärker skillnaden mellan spänningarna V+ och V- med en mycket hög förstärkning, A. Vi kan göra två antaganden om denna förstärkare. 1) Inresistansen är mycket hög vilket innebär att vi kan anta att strömmarna I+ och I- är lika med noll, så dessa behöver ni inte beräkna. 2) Utresistansen är mycket låg och

+

_

R1 R2 R3 R4

R5 R6

R7

I6

I7 I5

I4 I3 I2 I1

I+

I- V1

V2 V3 V4

V+

V-

Vut

(2)

detta innebär att spänningen Vut bestäms helt av skillnaden i inspänningar och inte alls av strömmen I7.

Nu kan vi gå vidare med att ställa upp ekvationer, och vi behöver 10 stycken eftersom vi har 10 stycken obekanta. Vi börjar från vänster i schemat och ställer upp följande ekvationer enligt Kirchoffs lagar och Ohms lag.

( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) ( 6 ) ( 7 ) ( 8 ) ( 9 ) ( 10 )

Nu är i princip problemet löst, det enda som är kvar är att implementera detta i ett Matlab program (skript). Börja med att skapa en m-fil i en texteditor (t.ex. medit eller emacs). Definiera sedan variabler för alla inparametrar och använd gärna kommentarer för att förklara vad som är vad. Använd följande värden för alla spänningar, resistanser och förstäkning:

% Matlab script för laborationsuppgift 1

% Inparametrar

% Spänningar V1 = 1; % V V2 = 2; % V V3 = 3; % V V4 = 4; % V

% Resistanser R1 = 1000; % Ohm R2 = 1000; % Ohm R3 = 1000; % Ohm R4 = 1000; % Ohm R5 = 2000; % Ohm

1 1 1

R V I =V − ⁺

2 2 2

R V I =V − ⁺

3 3

3 R

V I =V − ⁺

4 4 4

R V I =V − ⁺

5

5 R

I =V⁺

4 3 2 1

5 I I I I

I = + + +

6

6 R

I =V⁻

7

6 I

I =

7

7 R

V I =V^ut − ⁻

(

₊− ₋

)

= AV V V_ut

(3)

R6 = 1000; % Ohm R7 = 4000; % Ohm

% Förstärkning A = 100000;

Nästa steg är att skapa det ekvationssystem som skall räkna ut de obekanta variablerna. För att göra detta kommer vi använda Matlab:s matrismultiplikation, som fungerar på följande sätt för en 2x2 matris som multipliceras med en 2x1 matris:



 





⋅ +

⋅

⋅ +

= ⋅



 



⋅



 





21 22 11 21

21 12 11 11 21

11 22 21

12 11

b a b a

b a b a b

b a a

a a

Om vi nu skriver om samtliga ekvationerna på följande sätt:

( 11 ) ( 12 ) ( 13 ) ( 14 ) ( 15 ) ( 16 ) ( 17 ) ( 18 ) ( 19 ) ( 20 )

Notera hur alla termer med en obekant är placerade till vänster om likhetstecknet och övriga termer är placerade på andra sidan, som får värdet noll om inga andra termer finns. Om vi nu placerar alla obekanta i en vektor (motsvarar b i matrisexemplet) och faktorn framför den obekanta i en matris (motsvarar a i exemplet), så kan dessa 10 ekvationer skrivas i matrisform, se nedan.

1 1 1

1 R

V R I +V⁺ =

2 2 2

2 R

V R I +V⁺ =

3 3 3

3 R

V R I +V⁺ =

4 4 4

4 R

V R I +V⁺ =

0

5

5− ⁺ =

R I V

4 0

3 2 1

5−I −I −I −I = I

0

6

6 − ⁻ =

R I V

7 0

6 − I = I

0

7 7

7 − + ⁻ =

R V R I V^ut

=0 +

−AV₊ AV₋ V_ut

(4)













=













⋅













−

− +

0 0 0 0 0 0

1 0

0 0 0 0 0 0

1 1

0 1 0 0 0 0 0 0

0 0

0 1 1 0 0 0 0 0

0 1

0 0 1 0 0 0 0 0

0 0

0 0 0 1 1 1 1 1

0 0

1 0 0 1 0 0 0 0

0 0

1 0 0 0 1 0 0 0

0 0

1 0 0 0 0 1 0 0

0 0

1 0 0 0 0 0 1 0

0 0

1 0 0 0 0 0 0 1

4 4

3 3

2 2

1 1

7 6 5 4 3 2 1

7 7

6 5

4 3 2 1

R V

V V V I I I I I I I

A A

R R

R R R R

ut

Detta kan se lite mastigt ut till en början, men ta t.ex. första raden så ger matrismultiplikationen:

1 1 1 1 1

7 6 5 4 3

2

1 1 0 0

0 0

0 0 0 0

1 R

V R I V V V

R V I I I I

I I

I + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ _ut = + =

⋅ ₊ ₋ ⁺

Kontrollera gärna övriga rader, så får ni lite känsla för matrismultiplikationer.

Alla element i 10x10 matrisen och 10x1 matrisen till höger om likhetstecknet innehåller kända värden från de första raderna i m-filen. Den andra 10x1 matrisen innehåller endast obekanta och det är alltså dess element vi söker. En kompakt form av ekvationssystemet kan skrivas på följande sätt:

V X M⋅ =

Där M är 10x10 matrisen, X är 10x1 matrisen med de obekanta och V är den kända 10x1 matrisen. I er m-fil kan ni nu fortsätta med att lägga till rader som definierar de två matriserna M och V i Matlab, och kalla dem också för M och V. Exemplet nedan ger lite hjälp på traven. Var väldigt noga med hur storleken på era matriser blir, kolla gärna genom att köra ert program följt av kommandot whos.

% 10x10 Matris

M = [1 0 0 0 0 0 0 1/R1 0 0 ; ... % 1:a raden 0 1 0 0 0 0 0 1/R2 0 0 ; ... % 2:a raden 0 0 1 0 0 0 0 1/R3 0 0 ; ... % 3:e raden

% osv ....

När ni är klara med definitionen av M och V återstår bara att lösa ut vektorn med obekanta, som vi kallar för X. Detta gör ni genom att matrisdividera på båda sidor om likhetstecknet precis som det beskrevs i kurslitteraturen (Introduktion till Matlab) under sektionen Matematiska operationer. Var noggranna med att använda rätt typ av matrisdivision annars kommer det inte att fungera.

Uppgifter

Börja först med att beräkna lösningen för de inparametrar som var givna, och fyll i svaren från matrisen X i tabellen nedan.

Försök sedan lista ut hur utspänningen, Vut, beror av inspänningarna, V1, V2, V3 och V4, genom att testa er fram. Fyll i de tester ni gör i tabellen. Ett tips är att ändra en

(5)

inspänning i taget. När ni har kommit på hur formeln ser ut kan ni skriva den i fältet nedanför tabellen.

Redovisa också det program ni skrivit.

Tabell 1. Resultaten från beräkningarna. Observera att enheten på strömmarna skall anges i mA men i variabeln X har de enheten A.

V1

[V] V2

[V] V3

[V] V4

[V] I1

[mA] I2

[mA] I3

[mA] I4

[mA] I5

[mA] I6

[mA] I7

[mA] V+

[V] V-

[V] Vut

[V]

1 2 3 4

Skriv in formeln för Vut här:

(6)

Uppgift 2

Beskrivning

Att beräkna derivatan av en given funktion är en ganska vanlig uppgift vid numerisk analys av bl.a. mätdata. I denna uppgift ska vi göra en funktion som beräknar derivatan av en funktion då x-värdena och y-värdena är givna som vektorer.

En derivata kan beräknas numeriskt på lite olika sätt, den vanligaste metoden kallas för centrala differensmetoden och finns beskriven nedan.

x

x x f x x f dx x df

f ∆

∆

−

∆

= +

= 2

) (

) ) (

(

' där y= f(x)

Nu antar vi att vi har en viss mängd x-värden lagrade i en vektor (t.ex. en 20x1 matris) och vi kallar varje element i denna vektor för x1, x2, x3, ..., x20. Sedan beräknar vi lika många y-värden med funktionen f(x) och får då en 20x1 vektor med elementen y1=f(x1), y2=f(x2), y3=f(x3), ..., y20=f(x20). På samma sätt som formeln ovan kan vi med dessa vektorer beräkna den numeriska derivatan genom följande formel.

1 1

1

' 1

− +

−

= −

n n

n x x

y

y y där n=2,3,...,19

Denna formel har en fördel jämfört med den tidigare, eftersom den inte är beroende av att alla x-värden är jämnt fördelade. I den första formeln används ju ∆x som skillnad i x-värden. Observera att det uppstår lite problem för ändpunkterna med denna formel (n=1 och n=20), så där får vi använda följande formel.

1 2

'1

x x

y y y

−

= − och

19 20

'20

x x

y y y

−

= −

Men strikt numerisk ger denna formel derivatan mitt emellan x1 och x2 respektive x19

och x20, så ett litet fel får vi faktiskt här. Tyvärr går detta inte att göra så mycket åt, utan vi får leva med detta fel. I vissa speciella fall går det att få bättre noggrannhet, men då måste man ha någon kännedom om hur derivatan ser ut vid ändpunkterna.

Båda formler går att beräkna mycket snabbt och effektivt med Matlab genom att använda matrismultiplikation. Studera matrismultiplikationen nedan.













−

=













⋅













−

19 20

18 20

2 4

1 3

1 2

20 19 3 2 1

1 1 0 0 0

1 0 1 0 0

0 1

0 1 0

0 0

1 0 1

0 0

0 1 1

x x

x x x x x

M M

L L

O O O M

M L L

Resultatet blir faktiskt nämnarna i de båda formlerna ovan, och multiplicerar vi samma matris med y-värdena så får vi täljarna. Det enda som återstår är att utföra en elementvis division ( ./ )av de båda resultaten.

Skapa en m-fil som innehåller följande funktionshuvud.

% Funktion för att derivera y=f(x)

(7)

function dydx = derivera(x,y)

if size(x,2) > 1 x = x';

end

if size(y,2) > 1 y = y';

end

if size(x) ~= size(y)

error('vektorerna x och y måste ha samma storlek');

end

Spara sedan filen med namnet derivera.m. Detta funktionshuvud innehåller också en del felkontroll av argumenten. Första if satsen kontrollerar om vektorn x har fler än en kolumn, och i så fall byter den på kolumnerna och raderna genom att sätta tecknet ’ efter vektorn (kallas för att transponera en matris eller vektor). Nästa if sats gör samma sak med vektorn y. Den sista if satsen kontrollerar att x och y har samma storlek, om så inte är fallet visas ett felmeddelande. Kommandot size(A) ger alltså antalet rader och kolumner hos matrisen A, och om ett andra argument är givet för kommandot så anger den antalet element längs den dimension som anges (rader = 1 och kolumner = 2).

Det som är kvar för att ni skall kunna göra klart funktionen är att konstruera den matris som ger täljare och nämnaren, kalla t.ex. denna matris för D. Detta är ganska komplicerat att göra med hakparenteser ( [] ) som vi gjort tidigare eftersom storleken på vektorerna är okända just nu. En praktisk funktion att använda här är diag som skapar matriser från en vektor där vektorns element läggs på diagonalen i matrisen. Skriv ut hjälpen för detta kommando i Matlabfönstret.

>> help diag

DIAG Diagonal matrices and diagonals of a matrix.

DIAG(V,K) when V is a vector with N components is a square matrix of order N+ABS(K) with the elements of V on the K-th diagonal. K = 0 is the main

diagonal, K > 0 is above the main diagonal and K <

0 is below the main diagonal.

DIAG(V) is the same as DIAG(V,0) and puts V on the main diagonal.

DIAG(X,K) when X is a matrix is a column vector formed from the elements of the K-th diagonal of X.

DIAG(X) is the main diagonal of X. DIAG(DIAG(X)) is a diagonal matrix.

Example

(8)

m = 5;

diag(-m:m) + diag(ones(2*m,1),1) + diag(ones(2*m,1),-1)

produces a tridiagonal matrix of order 2*m+1.

Uppgift 3

På samma sätt som förra uppgiften, så går det också att beräkna integralen av en funktion given som x-värden och y-värden. Precis som för derivator så finns det flera sätt att beräkna integraler. Vi kan här bestämma oss för att använda trapetsmetoden, och den fungerar enligt formeln nedan, där vi antar att vi har en x-vektor och en y- vektor med 20 element vardera.

∫

¹ ⁼

1

0 ) (

x

dx x f

( )( )

( )

∫

² ⁼ ⁺ ⁻

1

1 2 2 1 21

) (

x

x x y y dx x f

( )( ) ( )( )

( )

∫

³ ⁼ ⁺ ⁻ ⁺ ⁺ ⁻

1

2 3 3 2 1 2 2 2 1

) 1

(

x

x x y y x x y y dx x f

osv. fram till sista elementen i x och y vektorerna.

( )( ) ( )( ) ( )( )

( )

∫

²⁰ ⁼ ⁺ ⁻ ⁺ ⁺ ⁻ ⁺ ⁺ ⁺ ⁻

1

19 20 20 19 2

3 3 2 1 2 2 2 1

) 1

(

x

x x y y x

x y y x x y y dx

x

f K

Dvs. integralen får också 20 element eftersom övre integrationsgränsen motsvarar varje element i x-vektorn. För att utföra detta med matrismultiplikation måste vi först dela upp multiplikationen av parenteserna.













−

⋅













+ +

+

+ +

+

⋅

19 20

2 3

1 2

20 19 3

2 2 1

3 2 2 1

2 1

0 0 0

0 0

0

2 1

x x

y y y

y y y y

y y

M L

O M

M

M L L

Kontrollera denna matrismultiplikation, så ni ser likheten med uttrycken ovan.

Vektorn i uttrycket kan ni få fram på liknande sätt som föregående uppgift, men med en lite annorlunda matris.













⋅













−

20 2 1

1 1 0

0

0 1

1 0

0 0

1 1

x x x

M L

O O M

M L

Notera att matrisen har en kolumn mer än antalet rader. Dessa matriser är lite mer komplicerade att skapa än helt fyrkantiga matriser (samma antal rader och kolumner).

Det enklaste är att först skapa en fyrkantig matris och sedan dra bort den del man inte vill ha. T.ex. kan matrisen ovan skapas med följande kommandon.

% Matrisen för att få dx

D = diag(-ones(1,size(x,1))) + diag(ones(1,size(x,1)- 1),1);

(10)

D = D(1:end-1,:);

Först skapas huvuddiagonalen fylld med –1:or, och sedan adderas den övre diagonalen fylld med 1:or. Sedan tas den understa raden bort med radangivelsen 1:end-1, dvs. alla rader från nummer 1 till sista minus ett. Alla kolumner väljs genom att ange ett kolon ( : ) i kolumnangivelsen.

Matrisen med elementen i y-vektorn kan skapas på följande sätt.













⋅













20 19 19 3 3 2 2 1

0 0

0 0 0

0 0

0 0 0

1 1 1

1

0 1

1 1

0 0

1 1

0 0

y y y y y y y y

L L

O O M M

M M L L

L O O M M

M L L

Även här måste ni trixa lite för att skapa matriserna. Matrisen med y-värdena är ej fyrkantig, så ni behöver använda samma teknik som för den förra formeln. För att få matrisen med ettor så kan ni använda ett kommando som heter tril, som hämtar den undre delen av en matris från diagonalen och ersätter övriga element med nollor.

Sedan kan ni manuellt sätta det översta elementet till noll, se kodsnutten nedan.

% Trapetsmatrisen för y T = tril(ones(size(x,1)));

T(1,1) = 0;

Här skapar kommandot ones(size(x,1)) en fyrkantig matris med endast ettor och tril ger endast den nedre delen och övriga element blir nollor.

Det gäller att ha koll på hur många rader och kolumner varje matris har, eftersom det inte helt tydligt framgår av exemplen ovan. Dock går det att lista ut från multiplikationerna. Kom ihåg regeln att resultatet av en matrismultiplikation har samma antal rader som den första faktorn och samma antal kolumner som den andra faktorn. Detta är en bra övning i att skapa och använda matriser på ett effektivt sätt i Matlab.

Försök nu själva skapa funktionen som integrerar vektorn y med avseende på vektorn x, och kalla funktionen för integrera.m. Utgå gärna från funktionen derivera.m och ändra där det behövs, t.ex. kan raderna för felkontrollen återanvändas. Det svåra är att översätta alla matrisoperationer ovan till Matlab kod, men tag hjälp av det tidigare exemplet. Titta också i kurskompendiet.

Uppgifter

Integrera funktionen y=sin(x) för x-värdena från -2π till 2π med avstånden 0.5 och 0.1. Jämför med den exakta integreringen, -cos(x), genom att plotta numerisk och exakt integrering i samma figur. Försök förklara eventuella skillnader (Tips: vad är det som saknas i den exakta integreringen?).

(11)

Uppgift 4

Beskrivning

I denna uppgift skall ni använda funktionen integrera.m för att beräkna volymen på en rotationssymmetrisk kropp, där radien bestäms av formeln r=1+₂¹sin(x) där x varierar mellan − och ^π₂ ³₂^π . Denna integral beräknas ju genom att ta arean för en cirkel och integrera denna längs x-axeln enligt formeln:

( )

∫ ∫

− −

+

2 =

3

2

2 3

2

2 2

2 1 1sin( )

π

πr dx x dx

I Matlab går det också smidigt att rita denna kropp med kommandot surf. Det åstadkommer ni genom att rita den yta som omger den kropp ni beräknat volymen på.

Vad ni kommer att behöva är då (x,y,z) koordinater för ett antal punkter längs denna yta, och för att få dessa måste ni utföra rotationen och beräkna dessa koordinater, se Fig. 2.

Fig. 2. Den rotationssymetriska volymen och beräkningen av dess y och z koordinater.

Skapa en ny m-fil och skriv alla kommandon i den eftersom det då blir lätt att ändra i kommandosekvensen. Låt x variera mellan − och ^π₂ ³₂^π med steg om 0.2, och låt vinkeln a variera mellan 0 och 2π med steg om 0.2. Skapa sedan en tillhörande matriser med kommandot meshgrid.

% Skapa matriser med meshgrid [X,A] = meshgrid(x,a);

Räkna sedan ut y och z koordinaterna från dessa matriser enligt formeln för radien ovan. Koden för detta finns nedan.

% Räkna ut Y och Z koordinaterna Y = (1+0.5*sin(X)).*sin(A);

Z = (1+0.5*sin(X)).*cos(A);

Z

Y

X

r r sin(a)

r cos(a)

Y Z

(12)

Observera användandet av punkt-gånger för att utföra allt element för element. Rita sedan hela kroppen med kommandot surf, vilket utför nästan samma sak som mesh men skapar en yta istället för ett rutnät.

% Rita volymen surf(X,Y,Z)

Kör er kommandofil och studera figuren, prova att rotera den genom att använda rotationsverktyget i figurfönstret (symbolen ). För att snygga till det hela så kan ni utföra följande kommandon.

% Skapa en ljuskälla l = light;

% Gör volymen lite snyggare

set(s,'edgecolor','none','facelighting','phong','facecolo r',[0 1.0 0]);

% Fixa till aspect ration, lägg till perspektiv och ta bort de vita fälten

axis('equal')

set(gca,'projection','perspective','visible','off');

Dessa kommandon och genom att vrida figuren lite, så kan ni få bilden på framsidan av kurskompendiet.

Till slut kan ni lägga till en titel på figuren med kommandot title, där texten skall innehålla volymen på kroppen. Anta att alla mått är angivna i cm, så får ni en enhet också. Volymen beräknar ni med den funktion ni skapade i förra uppgiften. En annan funktion ni kommer att behöva är num2str, som omvandlar en numerisk värde till en teckenmatris (sträng). Använd Matlab:s kommando help för att se hur ni använder dessa (Tips: för att få cm³ skrivet i titeln skriver ni ’cm^3’).

Uppgifter

Rita endast upp figuren i ett fönster och försök få den att se ut som den i kurskompendiet. Laborera gärna med de olika grafiska verktygen i figuren.