Formatera str¨ angar f¨ or utmatning - Formatkoder och kommandot sprintf 77

N¨ar vi vill skriva ut ett resultat fr˚an en m-fil till kommandof¨onstret s˚a kan vi anv¨anda kom-mandot num2str(x) som omvandlar en siffra till motsvarande textstr¨ang. Vi kan d˚a bygga upp en textstr¨ang som vi sedan skriver ut i kommandof¨onstret med kommandot disp:

out = [’Primfaktorerna ¨ar ’ num2str(f(1)) ’ och ’ num2str(f(2))] f¨oljt av disp(out);.

Problemet med den konstruktionen ¨ar att vi m˚aste veta fr˚an b¨orjan hur m˚anga tal vi vill skriva in i str¨angen. Ett mer flexibelt s¨att att g¨ora samma sak ¨ar att anv¨anda komman-dot sprintf. Kommankomman-dot sprintf l˚ater oss skriva in text och variabler till en textstr¨ang med anv¨andande av formateringsregler, vi kan sedan skriva ut str¨angen p˚a sk¨armen med anv¨andande av disp, precis som vanligt. Att beskriva hur formaterad utskrift g˚ar till kan bli v¨aldigt abstrakt, s˚a l˚at oss b¨orja med ett exempel och f¨ors¨oka f¨orst˚a det hela utifr˚an det exemplet. I kommandot

y=2;

x=3;

string = sprintf(’Svaret ¨ar %f eller kanske %g ’,y ,x) skapar sprintf en text-str¨ang som kallas string av den str¨ang som st˚ar till v¨anster om det f¨orsta kommatecknet i parantesen. Men varje g˚ang en formatkod - ordet som inleds med % - f¨orekommer s˚a ers¨attes formatkoden med en variabel h¨amtad efter det f¨orsta kommatecknet. Antag att y = 2 och x = 3, str¨angen string blir d˚a lika med ” Svaret ¨ar 2.000000 eller kanske 3”. Men varf¨or skrivs y ut med 6 decimaler och x utan decimaler? Det beror p˚a att formatkoden inte bara visar att man skall stoppa in v¨ardet av en variabel p˚a den platsen i str¨angen, utan ocks˚a hur den skall formateras. %f s¨ager att det tal som skall stoppas in ¨ar ett decimaltal, och om vi inte anger annat s˚a kommer det att skrivas med 6 decimaler. %g anger att talet skall skrivas ut p˚a ”l¨ampligt” s¨att - om det ¨ar ett tal st¨orre ¨an 10⁻⁴ och mindre ¨an 10⁶ skrivs det ut som ett decimaltal, men utan efterf¨oljande nollor. Ligger talet utanf¨or det intervallet skrivs det som ett decimaltal g˚anger en exponent av tio. Vill man alltid ha det senare formatet kan man anv¨anda formatkoden %e. Prova till exempel f¨oljande (f¨or att skriva ut ett procenttecken i en textstr¨ang skriver man %%):

y=2.17;

x=34567.98765;

string = sprintf(’x och y med formatet %%f %f %f ’,x,y) disp(string)

string = sprintf(’x och y med formatet %%g %g %g ’,x,y) disp(string)

string = sprintf(’x och y med formatet %%e %e %e ’,x,y) disp(string)

Vi kan p˚averka b˚ade hur stor bredd sifferf¨altet skall uppta i textstr¨angen och hur m˚anga decimaler talet skall skrivas ut med. Det g¨ors genom att skriva in siffror mellan procent-tecknet och formatkoden. Den f¨orsta siffran talar om vilken minsta bredd sifferf¨altet skall uppta, den andra siffran talar om hur m˚anga decimaler vi vill se utskrivna, mellan dessa siffror skrivs en punkt. S˚a betyder %7.3f att vi vill skriva ut ett decimaltal s˚a att totala bredden blir minst 7 tecken och att talet skall ha 3 decimaler. Prova

y=2.17;

x=34567.98765;

string = ...

sprintf(’x har v¨ardet %7.3f och y har %7.3f med formatet %%7.3f ’,x,y);

disp(string)

f¨or att se hur det fungerar. L¨agg m¨arke till dels att talet skrivs ut med st¨orre bredd ¨an den minimibredd vi anger om detta ¨ar n¨odv¨andigt f¨or att vi skall f˚a plats med hela talet inklusive det antal decimaler vi ber om, dels att om talet ¨ar kortare s˚a fylls det p˚a med blanktecken f¨or att minimikravet p˚a bredd skall uppfyllas (vill vi att siffrorna skall st˚a till v¨anster i f¨altet skriver vi %-7.3f ist¨allet). Oftast vill man ju inte ha dessa blanka tecken inuti str¨angen utan skriver %0.3f f¨or att ange antalet decimalsiffror, men ibland kan man vilja linjera upp siffror som skrivs ut i p˚a varandra f¨oljande rader, och d˚a kan en fix vidd vara anv¨andbar. Prova ocks˚a vad som h¨ander om ni anv¨ander formaten

string = sprintf(’x is%7.3g and y is %7.3g with format %%7.3g ’,x,y); och string = sprintf(’x is %7.3e and y is %7.3e with format %%7.3e ’,x,y);

Ut¨over e, f och g s˚a finns det ett antal andra formatkoder, prova help sprintf. Den som vi kan komma att beh¨ova i den h¨ar kursen ¨ar s som anger att vi vill stoppa in en hel str¨ang, oavsett hur l˚ang den ¨ar (den kan allts˚a ha l¨angden 1 n¨ar vi bara vill skriva ut en enskild bokstav).

Det som g¨or sprintf l¨amplig att anv¨anda f¨or att skriva ut resultat n¨ar vi inte vet hur m˚anga variabler vi vill skriva ut ¨ar att om det finns fler variabler till h¨oger om det f¨orsta kommat inom parantesen efter sprintf s˚a b¨orjar tolkningen om fr˚an b¨orjan. Fel anv¨ant blir det inte s¨arskilt snyggt, prova hur out=sprintf(’Resultatet ¨ar %g ’,x,y) ser ut. Ist¨allet f˚ar vi konstruera fler str¨angar, varav en bara best˚ar av det ok¨anda antalet variabler som vi vill skriva ut, f¨or att sedan l¨agga ihop str¨angarna till en l˚ang str¨ang som vi l¨agger ut med hj¨alp av disp:

x=1:7;

out1=sprintf(’Resultatet ¨ar’);

out2=sprintf(’ %g’,x);

out=[out1 out2 ];

disp(out)

Observera att str¨angen out2 inneh˚aller ett mellanslag. Vi ¨ar nu mogna att koda den sista delen av v˚art program, d¨ar vi skriver ut primfaktorerna:

out1=sprintf(’Talet %0.0f har %0.0f primfaktorer \n De ¨ar’, x, AntalFaktorer);

out2=sprintf(’ %g’,Faktor);

out=[out1 out2];

disp(out)

8.3 Avancerad grafik

I detta avsnitt samlar vi n˚agra kommandon som kan anv¨andas f¨or mer avancerad grafik, b˚ade 2D-grafik och 3D-grafik som vi kommer att beskriva i n¨asta avsnitt.

8.3. AVANCERAD GRAFIK 79 8.3.1 Fler grafikf¨onster - kommandot figure

Normalt ritas alla grafer i ett f¨onster som heter figure(1). Ibland n¨ar man till exempel analyserar resultaten fr˚an en laboration s˚a vill man k¨ora igenom ett program som g¨or hela analysen och spara ett antal olika plottar. Detta kan man g¨ora genom att ge kommandot figure(#) d¨ar # ¨ar ett tal. Alla plottar som ritas efter detta kommando kommer att ritas i ett grafikf¨onster som heter figure(#) ¨anda tills vi ger ett nytt kommando med ett nytt f¨onsternamn.

8.3.2 Att redigera en graf interaktivt

I ComsolScript kan vi redigera grafer f¨or hand direkt i det f¨onster d¨ar grafen ritas. Uppe i verktygslisten finns en knapp som ser ut som en liten palett. Trycker man p˚a den ¨oppnas ett f¨onster “Redigera graf”. F¨or att demonstrera hur man redigerar grafer k¨or vi f¨orst f¨oljande kod:

close all clear all

x=0:.2:4; y= sin(x); z=cos(x);

subplot(2,1,1) plot(x,y) subplot(2,1,2) plot(x,z,’r*’)

Innan du b¨orjar redigera m˚aste du v¨alja vilket objekt du vill p˚averka. Detta g¨or du genom att v¨alja n˚agot av de element som listas i v¨ansterspalten. I ComsolScript kallas det vi normalt kallar graf f¨or “Axes”. Eftersom vi har tv˚a grafer i grafikf¨onstret finns det tv˚a objekt som kallas “Axes” i v¨ansterspalten. Objekten “Axes” kommer i samma ordning som numreringen i subplot-kommandot, s˚a n¨ar vi klickar p˚a den f¨orsta “Axes” i listan, eller n˚agot objekt under det, s˚a p˚averkar vi enbart egenskaper i det f¨onster som svarar mot position ett osv.

V¨aljer vi “Axes” i v¨ansterspalten ¨oppnas ett f¨onster d¨ar vi kan p˚averka ett antal egenskaper hos grafen.

• Axis. H¨ar kan vi redigera axlarna och st¨odrastret. Normalt ¨ar “Axis” valt i den

¨

oversta flikmenyn n¨ar f¨onstret ¨oppnas. Vi kan d˚a v¨alja manuella v¨arden f¨or min och max p˚a skalan genom att vid ”x/y-limits” ta bort markeringen vid “Auto” och skriva in nya min och maxv¨arden. Vi kan ocks˚a v¨alja att ha logaritmisk skala p˚a en eller b¨agge axlar. Fyller vi i boxen “Axis equal” s˚a kommer l¨angdskalan p˚a x- och y-axeln att s¨attas lika.

• Title. H¨ar kan vi v¨alja en titel f¨or hela figuren, och etiketter f¨or x- och y-axeln.

• Font. H¨ar kan vi v¨alja typsnitt, storlek, stil och f¨arg f¨or texten i grafen.

Om vi ist¨allet v¨aljer “line” i v¨ansterspalten kan vi p˚averka hur data representeras i grafen.

Vi kan ¨andra f¨arg, linjebredd, linjetyp. Vi kan ocks˚a v¨alja vilken typ av mark¨or som visas i datapunkterna och skriva dit en dataetikett

8.3.3 Handtags-grafik

ComsolScript ger oss m¨ojlighet att i detalj via kommandon i en m-fil styra utseendet av den grafik vi producerar p˚a ungef¨ar samma s¨att som n¨ar vi redigerar interaktivt som ovan.

Detta kan vara av v¨arde n¨ar vi producerar grafik i ett program och vill p˚averka utseendet av graferna inifr˚an programmet ist¨allet f¨or att g¨ora det interaktivt. Detta g¨ors med hj¨alp av

n˚agot som kallas ”handtagsgrafik” (eng handle graphics). N¨ar vi ¨agnar oss ˚at datorprogram-mering s˚a ¨ar ett ”handtag” en variabel som pekar tillbaks p˚a ett objekt vars egenskaper vi vill kunna p˚averka. Det kan vara en fil vi har ¨oppnat i ett program och som vi vill skriva data till eller st¨anga eller byta namn p˚a, eller som i det h¨ar fallet ett diagram vars egenskaper vi vill kunna ¨andra. Handtagen skaffar man sig genom kommandot handtag = cga. Detta ger oss ett handtag som vi sedan kan anv¨anda dels f¨or att skaffa oss information om grafer, dels ¨andra variabelv¨arden, och d¨armed egenskaper hos grafen. Handtagsgrafik ¨ar i sig en sk¨on konst, vi kan inte g¨or mer ¨an skrapa p˚a ytan h¨ar, f¨or en mer fullst¨andig beskrivning h¨anvisar vi till ?/Plotting and Visualizing Data/Axes/Overview of Axes Functions Vi kan f˚a en lista ¨over de egenskaper vi f˚ar tillg˚ang till genom det handtag vi skapar n¨ar vi ger plotkommandot ¨ar linjens egenskaper, om vi skriver get(handtag) i kommandof¨onstret s˚a f˚ar vi en lista p˚a de egenskaper hos detta objekt som vi kan ¨andra. Tabell 9.7 i ?/Plot-ting and Visualizing Data/Axes/Axes properties ger en fullst¨andig beskrivning av vilka egenskaper hos grafikf¨onstret du kan ¨andra. Motsvarande egenskaper hos grafen sj¨alvt beskrivs i Tabell 9.11 i ?/Plotting and Visualizing Data/2D Graphics/Editing Plots L˚at oss genom ett exempel visa hur man kan arbeta med handtagsgrafik, b¨orja med att skapa f¨oljande plot:

clear all clf

x=0:.1:4; y=sin(x); z=cos(x);

subplot(2,1,1) plot(x,y) handtag1 = gca subplot(2,1,2) plot(x,z,’r-+’) handtag2 = gca

Prova nu att modifiera dessa plottar t.ex. p˚a f¨oljande s¨att:

C set(handtag1,’Ylim’,[-1. 1.2]) C xnu = get(handtag1,’Xlim’) C xnu(1)=xnu(1)-0.5

C set(handtag1,’Xlim’,xnu) C XTick=-0.5:.25:4.5;

C set(handtag1,’XTick’,XTick) C set(handtag1,’FontSize’,10)

C handtag2line = get(handtag2,’Children’)

% skaffar handtag till linjen i graf 2 C get(handtag2line)

% visar vilka objekt vi kan p˚averka

set(handtag2line,’LineStyle’,’–’,’LineWidth’,2,’LegendString’,’cosinus x’,’Legend’,’on’)

8.4. 3-D GRAFIK 81

8.4 3-D grafik

Tre-dimensionell grafik ¨ar ju n˚agot av en h¨agring, eftersom vi alltid f˚ar h˚alla tillgodo med tv˚ a-dimensionella projektioner av kroppar och ytor i den trea-dimensionella rymden. ComsolScript hj¨alper oss att skapa s˚adana h¨agringar p˚a flera olika s¨att. Vi skummar litet p˚a ytan, och h¨anvisar den som beh¨over en grundlig genomg˚ang till ?/Plotting and Visualizing Data/3D Graphics.

8.4.1 Kurvor i rymden - plot3

Kommandot plot3(x,y,z,s) fungerar ekvivalent med kommandot plot i tv˚a dimen-sioner. Om x, y och z ¨ar vektorer (med samma l¨angd) s˚a plottas kurvan mellan punkterna, str¨angen s styr linjens utseende och har samma format som f¨or plot-kommandot. Argumenten x, y och z kan ocks˚a vara matriser av samma storlek, kommandot plottar d˚a en kurva f¨or varje kolumn i matriserna. Vill man ha fler kurvor i samma graf, men med olika antal punkter, eller med olika egenskaper hos linjen kan man ge kommandot p˚a formen plot3(x1, y1, z1, s1, x2, y2, z3, s2,...). I ¨ovrigt fungerar de flesta kommandon som kan anv¨andas p˚a en tv˚adimensionell plot ocks˚a h¨ar:

clear all;

z=1:pi/32:2*pi;

y=sin(z);

x=cos(z);

plot3(x,y,z,’r’) xlabel(’cos(t)’) ylabel(’sin(t)’) zlabel(’t’) grid on

title(’Helix i tre dimensioner’)

I grafikf¨onstret kan vi rotera grafen f¨or att betrakta den ur valfritt perspektiv.

Tryck ned knappen “Kretsa/Panorera/Zooma”. Om vi trycker ned musknappen med mark¨oren n˚agonstans i grafikf¨onstret och drar kommer grafen att rotera. Pr¨ova!

Vi kan ocks˚a zooma in och ut med knappen “Glidning in/ut”. N¨ar man b¨orjar rotera grafen ¨ar det l¨att att g˚a vilse och till slut inte veta vad som ¨ar upp och ned. D˚a kan knappen “G˚a till standard 3D vy” hj¨alpa oss hem igen.

Beroende p˚a vilka inst¨allningar som g¨allde senast du var inne och plottade visas ibland den 3-dimensionella grafen utan axlar. Detta kan du ¨andra genom att ge box on. Alternativt kan du ge enbart box, som togglar, dvs s¨atter box=’on’ om den ¨ar ‘off’ och omv¨ant.

8.4.2 Funktionsytor - mesh och surf

Funktionsytor kan ˚ask˚adligg¨oras p˚a tv˚a s¨att, antingen som ett galler eller som en hel yta.

Innan vi ger oss h¨an m˚aste vi, i analogi med hur vi f¨orst skapade en x-vektor f¨or att g¨ora en tv˚adimensionell plot av en funktion skaffa oss ett x-y galler d¨ar funktionen skall utv¨arderas.

Detta g¨ors med kommandot [x, y] = meshgrid(xintervall, yintervall) som ger tv˚a matriser x och y som inneh˚aller koordinaterna f¨or varje punkt d¨ar funktionen skall utv¨arderas.

Prova till exempel:

clear all;

clf figure(1) xkoord=-8:.1:4;

ykoord=-4:.1:8;

[x, y]=meshgrid(xkoord, ykoord);

z=sin(sqrt(x.∧2+y.∧2));

mesh(x,y,z);

figure(2) surf(x,y,z);

8.4.3 Att styra utseendet av en 3D-graf Perspektiv

F¨orutom att l¨agga till axel-rubriker och titlar, och ¨andra axlarna gradering som f¨or en tv˚adimensionell graf (och naturligtvis hela batteriet med manipuleringar som handtags-grafiken ¨oppnar upp) s˚a finns det n˚agra speciella funktioner f¨or en 3D-graf. En ¨ar komman-dot view som anger ur vilken riktning vi betraktar en graf. Denna anges med tv˚a vinklar (i grader): theta och phi. Theta ¨ar vinkeln mellan x-axeln och det plan i vilket grafen lig-ger. Theta=0 anger att x-axeln l¨oper parallellt med bildytan, och att y-axeln pekar in˚at i grafen (ser inget vidare ut!). Phi-vinkeln anger den vinkel som synlinjen bildar med horison-talplanet, phi=0 anger att vi ser p˚a grafen horisontellt fr˚an sidan, phi=90 att vi betraktar det hela rakt ovanifr˚an. Dessa vinklar kan man ¨andra genom kommandot view(theta,phi).

Man kan ocks˚a ge view(‘xy’) osv om man vill visa plotten i xy-planet.

F¨argskala

F¨or att rita in en f¨argkarta som anger vilken f¨arg som svarar mot vilket v¨arde ger man kommandot colorbar. Tycker man inte att f¨args¨attningen duger s˚a finns det ett antal olika f¨argkartor att v¨alja p˚a - se tabell 9-10 i ?/Plotting and Visualizing Data/3D Graph-ics/Colormaps and Color Bars f¨or att se vilka m¨ojligheter som finns. F¨orutom att ange vilken typ av f¨argkarta du vill ha kan man ocks˚a ange mellan vilka max- och min-v¨arden f¨arskalan skall str¨acka sig. Du kan experimentera litet, som till exempel:

C surf(x,y,z,’ColorMap’,’bone(256)’,’Clim’,[-4 4]) C colorbar

C surf(x,y,z,’ColorMap’,’bone(256)’,’Clim’,[-1 1]) C colorbar

C surf(x,y,z,’ColorMap’,’jet(256)’,’Clim’,[-1 1]) C colorbar

C surf(x,y,z,’ColorMap’,’jet(256)’,’Clim’,[-2 2]) C colorbar

C surf(x,y,z,’ColorMap’,’jet(16)’,’Clim’,[-2 2]) C colorbar

f¨or att se hur det fungerar.

ColorMap och Clim ¨ar b¨agge egenskaper som kan ¨andras genom att anv¨anda handtagsgrafik.

En yta som plottas byggs internt upp av ett antal sm˚a element, pixlar. Hur dessa f¨argas styrs av kommandot shading. Det finns tre varianter, shading(‘flat’) visar varje en-skilt bildelement med en f¨arg, vald i en av h¨ornen av elementet. shading(‘faceted’¨ar samma sak, men med varje element omges av en svart kant. Vi kan slutligen f˚a en utj¨amnad f¨args¨attning utan rutn¨at med shading interp.

8.4. 3-D GRAFIK 83 Har vi gjort en snygg f¨argad yta kan det ibland vara effektfullt (men oftast inte s¨arskilt informativt) att helt fril¨agga den innan vi klipper in den i n˚agot annat program. Det g¨ors genom grid off och axis off.

8.4.4 Konturplottar och projektioner

Vi kan rita konturplottar med kommandot contour(x,y,z), antalet niv˚aer kan styras med ett fj¨arde argument.

8.5 Ovningsuppgifter ¨

1. Rita f¨oljande funktioner, prova med olika kombinationer av mesh, surf, surfl och shadint interp, samt testa olika f¨argkartor.

z = q

x²+ y² − 10 ≤ x ≤ 10, −10 ≤ y ≤ 10 z = sin

q

x²+ y²



− 2π ≤ x ≤ 2π, −2π ≤ y ≤ 2π

Kapitel 9

Ber¨ akningar och anpassningar

9.1 Enkel statistik

ComsolScript ¨ar i allm¨anhet mycket bra p˚a att hantera matriser och vektorer. Det l¨onar sig att f¨ors¨oka att konsekvent utnyttja det n¨ar vi behandlar data. Antag till exempel att vi gjort en serie m¨atningar (3.6, 5.4, 2.8, 8.3, 3.7, 4.2, 6.1, 4.1, 3.8, 2.3) och vill ber¨akna medelv¨arde och standardavvikelse f¨or dessa. Vi kan d˚a naturligtvis g¨ora detta genom att tilldela 10 olika variabler dessa v¨arden: x1 = 3.6, x2=5.4, x3=8.3 ... och sedan ber¨akna < x > = (x1 + x2 + x3 +...+x10)/10, och sedan forts¨atta p˚a motsvarande s¨att f¨or standardavvikelsen. Men det ¨ar mycket enklare att ist¨allet fylla en vektor med m¨atdata och sedan operera p˚a vektorn:

medel = sum(x)/length(x) medel =

4.4300

diff = x - medel

diff2 = diff.∧2

sigma = sqrt(sum(diff2)/(length(x)-1)) sigma =

1.7538

Annu enklare ¨¨ ar givetvis att anv¨anda de inbyggda funktionerna mean(x) och std(x), men innan ni g¨or det m˚aste ni en g˚ang ha kodat detta f¨or hand!

9.2 Polynom

ComsolScript har n˚agra inbyggda funktioner f¨or att hantera polynom, f¨or en fullst¨andig be-skrivning h¨anvisas du till ?/Data Analysis, Statistics, and I/O/ Interpolation and Polynomials/Working With Polynomials. Vi g˚ar h¨ar igenom n˚agra av de viktigaste funktionerna.

9.2.1 Hitta r¨otter till polynom.

Givet ett polynom som till exempel y = 3x³− 7x + 4, s˚a ber¨aknas r¨otterna (nollst¨allena) till polynomet med funktonen roots(p) d¨ar p ¨ar en vektor med koefficienterna f¨or varje term i fallande ordning. Notera att om n˚agon koefficient ¨ar noll s˚a m˚aste detta explicit anges med

85

en nolla p˚a r¨att st¨alle i vektorn p.

 p = [3 0 -7 4];

 roots(p) ans =

-1.7583 1.0000 0.7583

9.2.2 Finn polynomuttryck f¨or givna r¨otter

Om vi omv¨ant har ett antal nollst¨allen i en vektor n s˚a ger oss funktionen poly(n) det polynom som har dessa r¨otter:

 n =[-2 1 4 7];

 poly(n) ans =

1 -10 15 50 -56

vilket allts˚a uttyds som polynomet x⁴− 10x³+ 15x²+ 50x − 56.

9.2.3 V¨arden p˚a polynom

V¨ardet f¨or ett polynom med koefficientvektor p i punkten x ges av polyval(p,x):

m =[3 -2 4 1];

polyval(m,4.2) ans =

204.7840

Allts˚a f(4.2) = 204.7840 om f (x) = 3x³− 2x²+ 4x + 1 9.2.4 Derivator av polynom

Att derivera ett polynom ¨ar naturligtvis inte s¨arskilt sv˚art att g¨ora utan ComsolScripts hj¨alp.

Det kan dock ibland n¨ar man skriver generella program vara beh¨andigt att ha en funktion till hands f¨or detta, i ComsolScript heter den polyder(m).

 m = [3 -1 3]

 polyder(m) ans =

6 -1

S¨ager oss allts˚a att derivatan av 3x²− x + 3 ¨ar 6x − 1.

9.2.5 Produkter och kvoter av polynom

ComsolScript ber¨aknar produkter av tv˚a polynom med funktionen conv(a,b), samt kvoten mellan tv˚a polynom med funktionen deconv(a,b). F¨or mer information anv¨ander du hj¨alpfunktionen.

9.3 Matrisekvationer

Vi har m˚anga g˚anger i detta kompendium tjatat om hur bra ComsolScript ¨ar f¨or matris-ber¨akningar. Nu skall vi titta litet n¨armare p˚a detta f¨alt. Vi b¨orjar med n˚agra enkla exempel.

Inversen A⁻¹ till matrisen A ¨ar en matris s˚a att A⁻¹· A = 1. I linj¨ar algebra har vi l¨art oss

9.4. MINSTA KVADRATANPASSNING 87

metoder f¨or att ber¨akna inversen till matriser som A = 3 4

−2 1

!

. ComsolScript-anv¨andare kan sl¨oa litet och ist¨allet anv¨anda sig av kommandot inv:

A=[4 2; 1 0];

Eftersom den numeriska precisionen hos datorn ¨ar begr¨ansad till ett visst antal decimaler s˚a blir svaret oftast inte riktigt exakt 1 n¨ar man kontrollr¨aknar. M˚anga fysikaliska fenomen kan beskrivas av matematiska formler som l¨ampar sig f¨or att ber¨aknas i en matrisformalism genom att en matris som representerar en operation appliceras p˚a ett initialtillst˚and och resulterar i ett sluttillst˚and, vi kan skriva det som en matrismultiplikation: S = O · I. Vi kommer d˚a ofta att s¨oka v¨arden p˚a den operator O som resulterar i sluttillst˚andet S givet ett initialtillst˚and I. Det h¨ar kan l˚ata v¨aldigt abstrakt, och m˚anga g˚anger ¨ar det v¨al ocks˚a p˚a det viset. Men ber¨akningsm¨assigt handlar det om att vi k¨anner matrisen I och S (ofta kolumnvektorer) s˚a s¨oker vi matrisen O som satisfierar ekvationen ovan. Detta problem l¨oser vi genom att multiplicera med l¨ampliga inverser - fr˚an r¨att h˚all: S · I⁻¹ = O · I· I⁻¹ = O.

9.4 Minsta kvadratanpassning

Om vi antar att en upps¨attning m¨atpunkter (x_i, y_i) uppfyller sambandet y = a + bx s˚a kan vi med hj¨alp av maximum likelihood principen finna att den b¨asta (i meningen att summan av de kvadratiska avvikelserna minimeras) uppskattningen av parametrarna a och b ges av:

a =

Dessa ekvationer ¨ar relativt enkla att ber¨akna med hj¨alp av ComsolScript. L˚at oss exempli-fiera genom att l¨osa en variation p˚a uppgift 8.13 i Taylor, d¨ar vi skall anpassa m¨atdata givna i tabellen till uttrycket s = s₀ + vt.

Tid (s) -3 -1 1 3

Position (cm) 4.0 7.5 10.3 12.0 Os¨akerhet i 0.5 0.6 0.7 1.1 position (cm)

Felet i tiden ¨ar negligerbart. F¨or att utf¨ora anpassningen m˚aste vi ber¨akna ett antal ter-mer av typen ^Pwx vilket ¨ar relativt enkelt med ComsolScripts vektoroperationer och ele-mentvisa operatorer. Summan ¨ar ju uppbyggd s˚a att f¨orst skall varje element i w multipliceras med motsvarande element i x, d¨arefter skall dessa produkter summeras. Vektoroperationen wx = w .* x utf¨or den f¨orsta av dessa operationer och skapar en ny vektor d¨ar elementen

¨

ar produkten av motsvarande element i w och x vektorerna. Vi kan sedan enkelt summera dessa genom att applicera vektoroperatorn sum p˚a denna vektor. Vi f˚ar allts˚a n˚agot i stil

med

% linfit.m

% performs least squares fit to y=a+bx

% performs least squares fit to y=a+bx

In document ComsolScript 1.2 i kursen experimentella metoder (Page 83-0)