Numeriikkaa

MATLAB sopii hyvin erilaisten numeeristen ongelmien ratkaisuun. Alla muutamia perusongelmia ja niihin liittyviä MATLAB-esimerkkejä.

Yhtälöryhmien ratkaisu

Lineaarinen yhtälöryhmä ratkeaa MATLABin takakeno-operaatiolla. Jos yhtälöitä on enemmän kuin tuntemattomia (kerroinmatriisin aste on pienempi kuin tuntemattomien lukumäärä) saadaan pienimmän neliösumman ratkaisu.

Yhtälöryhmä

10x1-7x2	= 7
-3x1+2x2+6x3	= 4
5x1-x2+5x3	= 6

voidaan esittää matriisimuodossa Ax=b. MATLABissa muodostetaan ratkaisua varten kerroinmatriisi A ja oikean puolen vektori b:

» A=[10,-7,0;-3,2,6;5,-1,5]
A =
    10    -7     0
    -3     2     6
     5    -1     5
» b=[7,4,6]'
b =
     7
     4
     6

» x=A\b
x =
     0
    -1
     1

» A*x
ans =
     7
     4
     6

Interpolaatio ja datan sovitus

Tarkastellaan vektoreihin x ja y talletettuja datapisteitä. Polynomi-interpolaatio datapisteille saadaan polyfit-funktiolla. Se palauttaa sovituspolynomin kertoimet, joita vastaavan polynomin arvoja voi laskea polyval:n avulla.

Interpoloituja taulukkohakuja voi tehdä interp1- ja interp2-funktioilla.

Esimerkki: Valitaan funktiosta

tasavälisesti pisteitä ja piirretään polynomi-interpolaatiot, jotka kulkee valittujen pisteiden kautta. (Esimerkki on peräisin Runge-nimiseltä matemaatikolta (1901). )

rungefun=@(x)1./(1+25*x.^2)
x=linspace(-1,1);
plot(x,rungefun(x),'k');
hold on; axis([-1,1,-0.4,1.3])

for n = [4,6,8,10,20]
  xp = linspace(-1,1,n);
  yp = rungefun(xp);
  plot(x,polyval(polyfit(xp,yp,n-1),x));
end

hold off; axis normal;

Datan sovitus annettuun funktioluokkaan (malliin) on sukua interpoloinnille. Polynomin voi sovittaa, joko tarkasti tai pienimmän neliösumman mielessä.

Esimerkki: suoran y = a + b x sovitus datapisteisiin. Generoidaan pisteet x = 0,1,2,...10 ja y = 1 + 2x + epsilon, missä epsilon on nomaalijakautunut satunnaismuuttuja, odotusarvona 0 ja hajontana 0.5. Sovitetaan pienimmän neliösumman regressiosuora tähän pistejoukkoon.

x = linspace(0,1,10)'; 
y=2*x+1 + randn(10,1)*0.5;
A = [ones(length(x),1) x]
b=A\y
xx=linspace(min(x),max(x));
plot(x,y,'o',xx,b(1)+xx*b(2))

Splini-interpolaatiossa sovitetaan datapisteiden kautta kulkemaan paloittain määritelty käyrä, joka on kahden vierekkäinen x-datapisteen välissä kolmannen asteen polynomi ja jossa lisäksi asetetaan ensimmäiset ja toiset derivaatat datapisteissä samoiksi. Näin saadaan sovitus, joka usein näyttää kulkevan jouhevasti datapisteiden kautta. Splinisovitus saadaan komennolla spline, edellisen esimerkin datapisteisiin yksinkertaisimmillaan komennolla

» plot(x,y,'o',xx,spline(x,y,xx))

Epälineaariset yhtälöt

Epälineaariset yhtälöt ja yhtälöryhmät eivät ratkea yksinkertaisilla matriisioperaatioilla vaan tarvitaan iteratiivisia ratkaisijoita. MATLABissa reaalifunktion nollakohtaa voi etsiä funktion fzero avulla. Polynomiyhtälöä varten on funktio roots.

Piirretään funktio 9sin(x)-x ja valitaan hiirellä klikkaamalla alkuarvaus nollakohdalle. Kootaan alla olevat komennot skriptiksi plotscript1.m

%% plotscript1.m
clf                         % Clear Graphics
f=@(x) 9*sin(x)-x; % Funktiomääritys 
fplot(f,[0 10]); grid on;shg
title('Valitse hiirellä piste nollakohdan läheltä')
[x0,y] = ginput(1);
hold on
plot(x0,0,'o');
x00=fzero(f,x0);
display(['Nollakohta = ' num2str(x00)])
plot(x00,f(x00),'r*')    % r='red'
legend('Kuvaaja','Alkupiste x0','Nollakohta x00')

>> plotscript1
Nollakohta = 2.8226  % Tulos

Numeerinen integrointi

Reaaliarvoisen funktion f(x) määrätyn integraalin

Sinifunktion integraali välillä [0, /2]

» integral(@sin,0,pi/2)

ans =

    1.0000

Integroidaan funktio x cos(30x) sin(x) välillä [0,2]. Tehdään MATLAB-funktion sisään globaali laskuri, jotta nähdään kuinka monta kertaan funktiota on kutsuttu. Lisäksi talletetaan laskentapisteet ja piirretään kuva.

Tiedosto: fun.m
-------------------------
function y=fun(x)
global samples count
samples= [samples,x];
count=count+length(x);
y=x.*cos(30*x).*sin(x);

» global count samples
» count=0; samples=[];
» I=integral(@fun,0,2*pi)

I =

    0.0070

» count

count =

        1590

» fplot(@fun,[0,2*pi])
» hold on; plot(samples,fun(samples),'or'); hold off

Differentiaaliyhtälöt

MATLABin perusasennuksessa on mukana laskentarutiinit tavallisten 1. kertaluvun differentiaaliyhtälöryhmien alkuarvotehtävien ratkaisemiseksi. Fuktiot ode45, ode23 ja ode15s. Korkeampaa kertalukua olevat tavalliset dy:t on muutettavan 1. kl:n muotoon ennen ratkaisua. Osittaisdifferentiaalityhtälöiden ratkaisijoita on Differential Equation Toolbox:ssa

Esimerkki: Toisen kertaluvun tehtävä putoavasta kappaleesta, jossa ilmanvastus on suoraan verrannollinen nopeuteen.

MATLAB-ratkaisua varten kirjoitetaan oma, diffyhtälö(systeemi)n y'=f(t,y) määrittelevä (vektori)funktio.

function ydot=omafun(t,y)

missä t on riippumaton muuttuja (usein aika) ja y riippuva (vektori)muuttuja. Funktio palauttaa arvon/arvot f(t,y).

function yp=pallo_putoaa(t,y)
m = 0.125; g=9.81;
yp(1) = y(2);
yp(2) = -g-y(2)/m; 
yp=yp(:); % yp:n oltava pystyvektori

MATLAB-ode-ratkaisijan syntaksi on

» aika = linspace(0,4);
» H = 10; v0=0;
» [t,y]=ode45(@pallo_putoaa,aika,[H,v0]);
» plot(t,y(:,1))

Optimointi

MATLABin optimointirutiineita ovat fminsearch ja fminbnd. Erillisessä Optimization Toolbox:ssa on lisää rutiineita.