Dibuixar una funció
============

En aquesta pràctica veurem com es pot representar gràficament una funció (d'una variable) usant Matlab.

Recordem:

* Exemple:

* Exemple: Pintar la gràfica de la funció sin(x) en [-10,10] amb pas 0.2

x=-10:0.2:10; % x es un vector de 101 punts en [-10,10]
y=sin(x); % y es el vector de les imatges de x
plot(x,y); % crea una finestra on es pinta la gràfica de la funció via la taula (x,y)
grid; %(opcional) afegeix una graella al dibuix

- Comentari: La gràfica que pinta plot(x,y) és la taula de punts unida per trossos de rectes.

* Exemple: Us modificadors de la funció plot per la mateixa taula de punts de la funció sinus que abans. Per més detalls sobre possibles modificadors vegeu el document Funcions.html del numfactory.

plot(x,y,'ro'); % r = red; o = circle
hold on; % perqè no es tanqui la finestra i poder dibuixar més coses sobre el dibuix anterior
plot(x,y,'b-.'); % b = blue; -. = punts i ratllles
hold off; % quan fem el proper plot es tancarà la finestra i se'n obre una de nova

Usar l'intrucció inline és la forma més fàcil de definir noves funcions en matlab i usar-les.

* Exemple: Introduim la funcio f(x)=10*(1-exp(-x/3)*sin(10*x)) com a funció inline i la dibuixem entre 0 i 12 amb pas 0.01

fun=@(t) 10*(1-exp(-t/3).*sin(10*t)); % funció inline, t és variable de la funció
x=0:0.01:12; % x vector on la 1a component és 0, la darrera 12 i estan separades de 0.01 en 0.01
y=fun(x); % y valors de la funció fun en la taula x
plot(x,y);

- Comentari: Podem triar els noms que volguem per la funció i la variable. Per exemple, haguessim pogut fer:

g=@(x) 10*(1-exp(-x/3).*sin(10*x));
x=0:0.01:12;
y=g(x);
plot(x,y);

- Comentari (FONAMENTAL): El punt que hem possat en definir la funció inline, exp(-t/3).*sin(10*t), és bàsic en Matlab! Si no el posem els càlculs no van bé.

* Exemple: si volem definit inline la funcio f(x)=exp(1/x)*sin(7*x)/(1+x^2) cal fer:

f=@(x) exp(1./x).*sin(7*x)./(1+x.^2);

La funció find permet conèixer quines components d'un vector (o matriu) compleixen una certa condició. La sintaxi és:

indexos_del_vector = find( condició )

(per més informació feu: doc find)

Exemples d'us de la instrucció find aplicada a les taules de valors x i y que hem generat amb la funció fun definida inline. Recordem: y=fun(x) on x es la partició de punts de [0:12] triada.

* Exemple 1: Volem saber quins són els punts de la taula (x,y) pels quals la funció fun queda per sobre de 3 i quants són en total.

index = find(y > 3);
numTot = size(index,2)

Concretament: via la comanda find obtenim un vector index que ens diu quins són els punts de la taula y pels quals el seu valor és per sobre de 3.
P. ex., si treiem el ; de index = find(y > 3) veiem que les primeres component del vector index són:

1, 2, 3, 4, 5, 6, 7, 8, 9, 24, 25, etc.

Això vol dir que els següents valor de la taula y són per sobre de 3:

y(1), y(2), y(3), y(4), y(5), y(6), y(7), y(8), y(9), y(24), y(25), etc

i el reste (no llistats) per sota o igual a 3.
(Recordeu que p.ex., y(1) és el valor de fun en x(1): y(1)=fun(x(1))).
Aleshores, el que tenim és que les 11 primeres components del vector index són:

index(1)=1, index(2)=2, index(3)=3, index(4)=4, index(5)=5, index(6)=6, index(7)=7, index(8)=8, index(9)=9, index(10)=24, index(11)=25.

El que hem fet per obtenir quantes components té el vector index és el següent. Si d'executa:

size(index)

s'obté un vector de dues components que ens dóna el nombre de files i columnes de index (entés com una matriu). Com, de fet, index és un vector filera, la primera component de size(index) és 1 i la segona és el nombre de columnes/components del vector index:

size(index,1)
size(index,2)

* Exemple 2: Com que en aquest cas els valors de la taula y no són nombre enters (sinó float o double) la condició que posem dins de find no pot ser una igualtat (sí que val usar la igualtat en find si la taula de punts de y fos a valors enters).
Raó: estem avaluant una funció fun (de variable contínua) en un conjunt discret de punts i per tant fora miraculos que la igualtat triada és donés per algún d'aquests valors concrets de la discretització.
Per exemple pel dibuix veiem que la funció fun talla dos cops l'altura 18 en l'interval [0.12], però per cap punt de la taula de valors val exactament 18. Així, si fem (ull: la igualtat és demana fent y==18 i no pas y=18)

index = find(y == 18);
numIgual18 = size(index,2)

obtenim numIgual18=0.

* Exemple 3: Si ara fem que la condició de find sigui:

index = find(abs(y-18) < 0.15);
numIgual18 = size(index,2)

Obtenim que hi ha 2 punts del vector y pels quals la seva distància a 18 és mes petita que 0.15. Això vol dir que hi ha dos valors de x de la taula de punts triada en l'interval [0,12] pels quals el valor de fun en aquests punts dista de 18 menys que 0.15.

* Exemple 4: Nombre de punts del vector y que són més grans o iguals que 18.

index = find(y >= 18);
numMesGran18 = size(index,2)

Animació: Exemple de simulació d'un oscil·loscopi
===============================

Simularem el comportament d'un oscil·loscopi a partir de la gràfica de la funció sin(x) desplaçarem al llarg de l'eix x mitjançant una fase F (àngle) que anirem incrementant.
En l'animació escriurem per pantalla el valor de F que correspon a cada frame i pintarem els eixos de coordenades.
Concretament pintarem successivament la funcio sin(x+F) per diferents valors de la fase F entre -endFase i endFase que anem incrementant amb pas incFase.

endFase=7; % L'interval de variació de la fase F és [-endFase,endFase].
incFase=0.05; % Increment de la Fase
x=-2*pi:0.2:2*pi; % Taula de valors de x que pintarem en cada dibuix (frame)

% Per pintar els diferents frames farem un bucle for que repetira la matrixa instrucció per diferents fases F entre els valors -endFase i endFase, incrementats de incFase en incFase.

for F=-endFase:incFase:endFase % inici bucle for: fem variar al fase F en el rang triat.
y=sin(x+F); % No cal introduir sin inline. El Matlab ja la té.
plot(x,y,'r')
% ----------------------(El següent bloc és opcional)
axis([-2*pi 2*pi -1.2 1.2]);
grid;
line([-2*pi 2*pi],[0 0],'Color','k','LineWidth',2);
line([0 0],[-1.2 1.2],'Color','k','LineWidth',2);
text(1.15,1.1,['Fase = ' num2str(F,'%.4f')]); % escriu 4 decimals Fase
% ----------------------(Fi bloc opcional)
drawnow; % Repeteix un i altre cop el dibuix.
%pause(0.01) % opcional: per retardar l'execució (en segons)
end % fi bucle for

Podem guardar un video de l'animació anterior, p. ex. en format avi, afegint 5 linies a l'anterior (tot seguit el reproduim sencer i el salvem):

endFase=7;
incFase=0.05;
x=-2*pi:0.2:2*pi;

% obrim el fitxer de video:

v = VideoWriter('sinusoide.avi');
open(v);

% repetim els càlculs que ara es guarden en l'arxiu triat en el format triat.

for F=-endFase:incFase:endFase
y=sin(x+F);
plot(x,y,'r')
axis([-2*pi 2*pi -1.2 1.2]);
grid;
line([-2*pi 2*pi],[0 0],'Color','k','LineWidth',2);
line([0 0],[-1.2 1.2],'Color','k','LineWidth',2);
text(1.15,1.1,['Fase = ' num2str(F,'%.4f')]); % text(x,y,missatge)
drawnow;

% guardem cada frame de video:
writeVideo(v,getframe);

end

% tanquem el fitxer de video:
close(v)

* Exemple addicional: Pintar les òrbites de Lissajou:

x=A*cos(a*t+delta),
y=B*sin(b*t)

on t és el temps i A,a,delta,B,b valors donats de forma que a/b és racional.

Triem per exemple: delta=pi/2, A=1 i B=1 fixos, b=a+1, i pintem l'animació que s'obté variant a de 1 a 100 amb pas 1 (si el pas per a és 1, és el mateix escriure a=1:1:100 o a=1:100, ometent el 1 del mig pel pas que ja se sobre-enten que és 1).

delta=pi/2.; A=1; B=1;
for a=1:100
b=a+1;
t=0:0.1/b:2*pi; % instants de temps que pintem (diferents per cada b)
% Ara les taules x,y que pintem es generen les dues en termes de la taula t.
x=A*sin(a*t+delta);
y=B*sin(b*t);
plot(x,y);
drawnow;
pause(0.1)
end