Rellenar un área de color

Utilizamos los vectores en combinación con los operadores lógicos para producir ciertos efectos que nos pueden ser útiles en las representaciones gráficas.

El comando fill rellena un recinto cerrado del color especificado. El recinto está descrito por dos vectores xx e yy que contienen las abscisas y ordendas de los puntos del recinto. Estos vectores están formados por elementos y por porciones de otros vectores extraídos mediante operadores lógicos.

Primero, vamos a recordar como se extrae un vector de otro mediante operadores lógicos

>> x=linspace(-1.5,4.5,10)
x =-1.5000 -0.8333 -0.1667 0.5000 1.1667 1.8333 2.5000 3.1667 3.8333 4.5000
>> u=x>1 & x<3
u =     0     0     0     0     1     1     1     0     0     0
>> x(u)
ans =    1.1667    1.8333    2.5000

Como vemos, se extraen aquellos elementos cuyo índice se corresponde con el valor uno.

Dibujamos la función -x²+3x+4 y definimos la región comprendida entre la curva, el eje X y las rectas x=1 y x=3 que queremos colorear con el comando fill. La definición de la región es:

Punto (1,0), seguido de (1, f(1)), todos los puntos de abscisa 1<x<3 y sus correspondientes ordenadas, el punto (3, f(3)) y finalmente, el punto (3,0). En el vector xx guardamos las abscisas y en el vector yy, las ordenadas. Llamamos al comando fill y le pasamos los dos vectores y el color de relleno.

x=linspace(-1,4,50);
f=@(x) -x.^2+3*x+4;
y=f(x);
hold on
plot(x,y,'b')
line([-2 5],[0 0],'color','k'); %eje horizontal X
xx=[1 1 x(x>1 & x<3) 3 3];
yy=[0 f(1) y(x>1 & x<3) f(3) 0];
fill(xx,yy,'y'); %rellena un área de color especificado
title('área')
xlabel('X')
ylabel('Y')
hold off

Cálculo de impuestos

La fórmula para calcular los impuestos de los salarios brutos de los empleados de una hopotética empresa es la siguiente:

Calcular e imprimir en una tabla el salario real de los empleados cuyo salario bruto es: 800, 2300, 3000, 3800 y 4200.

sb=[800 2300 3000 3800 4200]; %salario bruto
for k=sb
    if k<1000
        descuento=0.05*k;
    elseif k<2500
        descuento=50+0.15*(k-1000);
    elseif k<4000
        descuento=275+0.3*(k-2500);
    else
       descuento=725+0.4*(k-4000); 
    end
    disp([k k-descuento])
end

Corremos el script en la ventana de comandos

>> salario
   	800         760
        2300        2055
        3000        2575
        3800        3135
        4200        3395

Podemos escribir este script de forma alternativa, sin necesidad de sentencias condicionales (if, elseif) e iterativas (for).

sb=[800 2300 3000 3800 4200]; %salario bruto
descuento=(sb<1000).*sb*0.05+...
    (sb>=1000 & sb<2500).*(50+(sb-1000)*0.15)+...
    (sb>=2500 & sb<4000).*(275+(sb-2500)*0.3)+...
    (sb>=4000).*(725+(sb-4000)*0.4);

    disp([sb' (sb-descuento)'])

Como podemos apreciar cada línea que calcula el vector descuento, es el producto elemento a elemento de dos vectores. El de la izquierda es un vector lógico cuyos elementos son ceros y unos. En la siguiente tabla se muestra su significado

Ahora veamos el producto elemento a elemento de los vectores (sb<1000) y sb.

>> (sb<1000).*sb
ans =   800     0     0     0     0

multiplicando el vector [800,0,0,0,0] por 0.05 obtenemos el vector [40,0,0,0,0] que corresponde al descuento para los salarios brutos menores de 1000 euros. Del mismo modo, podemos probar la siguiente línea que calcula el descuento de los salarios comprendidos entre 1000 y 2500 euros. El salario bruto en este intervalo es 2300 euros, por tanto su descuento es 50+(2300-1000)*0.15=245.

Los vectores de la izquierda y de la derecha son, respectivamente

>> (sb>=1000 & sb<2500)
ans =     0     1     0     0     0
>> (50+(sb-1000)*0.15)
ans =    20   245   350   470   530

El producto elemento a elemento de ambos vectores es [0,245,0,0,0] que corresponde al descuento del segundo salario bruto. Así, para cada vector de salarios brutos sb obtenemos el correspondiente vector de descuentos debidos a los impuestos.

Al correr el script obtenemos una tabla cuya primera columna son los salarios brutos y la segunda los mismos con la deducción por impuestos

        800         760
        2300        2055
        3000        2575
        3800        3135
        4200        3395

Binomio de Newton

$$\begin{gathered} {\displaystyle {(x+a)}^n=\left(\begin{array}{l}n\\ 0\end{array}\right)x^n+\left(\begin{array}{l}n\\ 1\end{array}\right)x^{n-1}a+\left(\begin{array}{l}n\\ 2\end{array}\right)x^{n-2}{a}^2+\mathrm{...} \\ +\left(\begin{array}{l}n\\ n-2\end{array}\right)x^2{a}^{n-2}+\left(\begin{array}{l}n\\ n-1\end{array}\right)x·a^{n-1}+\left(\begin{array}{l}n\\ n\end{array}\right)a^n} \end{gathered}$$

Definir una función denominada pol_newton que devuelva un vector que contenga los coeficientes del polinomio resultado del desarrollo del binomio, cuando se le pasa el grado n y el valor del término a.

Primero calculamos el número combinatorio ⁿC_m o n sobre m.

${\displaystyle \begin{array}{l}\left(\begin{array}{c}n\\ m\end{array}\right)=\frac{n!}{m!(n-m)!}=\frac{n·(n-1)·(n-2)····(n-m+1)}{1·2····(m-1)·m}\\ \left(\begin{array}{c}n\\ m\end{array}\right)=\left(\frac{n}{1}\right)\left(\frac{n-1}{2}\right)····\left(\frac{n-m+2}{m-1}\right)\left(\frac{n-m+1}{m}\right)\end{array}}$

Definimos la función combinatorio

function r=combinatorio(n,m)
    k=1:m;
    t=(n-k+1)./k;
    r=prod(t);
end

Alternativamente, podemos utilizar la función de MATLAB nchoosek tal como se muestra en la porción de código

>> combinatorio(4,2)
ans =     6
>> nchoosek(4,2)
ans =     6

Definimos la función pol_newton, pasándole la potencia n y el valor del coeficiente a, nos devuelve los coeficientes p del polinomio.

function p=pol_newton(n,a)
    p=zeros(1, n+1);   %reserva espacio en memoria e inicializa a cero
    for k=0:n
        p(k+1)=combinatorio(n,k)*a^k;
    end
end

Probamos la función pol_newton en la ventana de comandos

>>pol_newton(2,1)
ans = 1     2     1
>> pol_newton(3,1)
ans = 1     3     3     1
>> pol_newton(4,1)
ans = 1     4     6     4     1
>> pol_newton(3,-1)
ans =     1    -3     3    -1

Sumas y productos

Dado un conjunto de datos x₁,x₂....x_n, calcular la suma de los siguientes productos

s=(x₁-x₂)·(x₁-x₃)...(x₁-x_n)+(x₂-x₁)·(x₂-x₃)...·(x₂-x_n)+....+(x_n-x₁)·(x_n-x₂)...·(x_n-x_n-1)

En el capítulo Vectores aprendimos a extraer elementos de un vector con otro vector

>> x=[0.97 1.12 2.92 3.00];
>> indice=[1 3 4];
>> x(indice)
ans =    0.9700    2.9200    3.0000

Vemos se obtiene el vector x sin su segundo elemento

Cuando x tiene muchos elementos podemos extraer todos los elementos de x excepto el que está en una posición k determinada del siguiente modo.

>> x=[0.97 1.12 2.92 3.00 3.33 3.97 6.10 8.39 8.56 9.44];
>> n=length(x);
>> k=3;
>> indice=[1:k-1 k+1:n]
indice =     1     2     4     5     6     7     8     9    10
>> x(indice)
ans =  0.9700  1.1200  3.0000  3.3300  3.9700  6.1000  8.3900  8.5600  9.4400
>> k=1;
>> indice=[1:k-1 k+1:10]
indice =     2     3     4     5     6     7     8     9    10
>> x(indice)
ans =  1.1200  2.9200  3.0000  3.3300  3.9700  6.1000  8.3900  8.5600  9.4400

Ahora, podemos codificar fácilmente el resultado de la suma de los productos anteriores, recordando que la función prod obtiene el producto de todos los elementos de un vector.

x=[0.97 1.12 2.92 3.00 3.33 3.97 6.10 8.39 8.56 9.44];
suma=0;
n=length(x);
for i=1:n 
       suma=suma+prod(x(i)-x([1:i-1 i+1:n]));
end
fprintf('La suma es %2.3f\n ',suma)

Corremos el script en la ventana de comandos

>> prueba
La suma es 274113.275

Posición del asteriode Pallas

Gauss estaba interesado en calcular de un modo preciso la órbita del asteroide Pallas, a partir de las observaciones de su posición.

Disponía de 12 pares de datos (t_i,x_i) que deseaba interpolar con la suma de seis funciones armónicas.

x= A₀+A₁cos(ωt)+ B₁sin(ωt)+ A₂cos(2ωt)+ B₂sin(2ωt)+...+ A_mcos(mωt)+ B_msin(mωt)

Donde los coeficientes A₀, A₁...A_m, B₁...B_m valen, véase al final de la página Ajuste de datos con funciones armónicas.

${\displaystyle \begin{array}{l}{A}_0=\frac{{\displaystyle \sum _{i=1}^N{x}_i}}{N}\\ A_k=\frac{2}{N}{\displaystyle \sum _{i=1}^N{x}_i\cos (k\omega t_i)}{B}_k=\frac{2}{N}{\displaystyle \sum _{i=1}^N{x}_i\sin (k\omega t_i)}k=\mathrm{1,2,...}m\end{array}}$

En primer lugar, representamos los datos

t=0:30:330; 
x=[408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
hold on
plot(t,x,'ro','markersize',2,'markerfacecolor','r')
xlim([0 360])
xlabel('Ascensión (Grados)')
ylabel('Declinación (Minutos)')
title('Posición del asteroide Pallas')
grid on

Calculamos ahora los coeficientes A₀, A₁...A₆, B₁...B₆ mediante las expresiones anteriores. La primera intención es utilizar un doble bucle: uno para obtener las sumas parciales y otro para obtener cada uno de los seis coeficientes a y b.

n=length(x);
w=2*pi/360;
a0=sum(x)/n; 
a=zeros(1,6); %reserva espacio en memoria para seis elementos
b=zeros(1,6);
for k=1:6
    for i=1:n
       a(k)=a(k)+x(i)*cos(k*w*t(i))*2/n;
       b(k)=b(k)+x(i)*sin(k*w*t(i))*2/n;
   end
end

Si tenemos en cuenta como se realizan las operaciones con matrices, podemos reducir significativamente las líneas de código marcadas en negrita:

• El producto de dos matrices de dimensiones (m×k)×(k×n)=(m×n) es una matriz de dimensión m×n
• El producto de un escalar por una matriz es otra matriz de la misma dimensión, cuyos elementos han sido multiplicados por el escalar.
• El coseno o el seno de una matriz, por ejemplo cos(A) es la matriz B(i,j)=cos(A(i,j)), se calcula el coseno de cada uno de los elementos de la matriz, obteniéndose otra matriz de la misma dimensión.

k=1:6;
a=cos(w*k'*t)*x'*2/n
b=sin(w*k'*t)*x'*2/n

k es un vector fila de seis elementos, es decir una matriz de dimensión 1×6. Su traspuesta k' es una matriz de dimensión 6×1. t es un vector fila de 10 elementos, de dimensión 1×10. El producto k'*t es una matriz de dimensión 6×10. Al calcular el coseno, se calcula el coseno de cada unos de los elementos de dicha matriz, la dimensión de la matriz resultante 6×10 no cambia .

x es un vector fila de dimensión 1×10, por tanto el producto con su traspuesta x', cos(w*k'*t)*x' es una matriz de dimensión (6×10)×(10×1)=6×1.

Al correr el script, obtenemos dos vectores columna a y b de seis elementos.

a =
 -411.0144
   43.4167
   -4.3333
   -1.0833
    0.3477
    0.1667
b =
 -720.2279
   -2.1651
    5.5000
   -1.0104
   -0.2721
   -0.0000

Conocidos los coeficientes A₀, A₁...A₆, B₁...B₆, representamos la función x(t), junto a las medidas (t_i,x_i) en la misma ventana gráfica (hold on)

x= A₀+A₁cos(ωt)+ B₁sin(ωt)+ A₂cos(2ωt)+ B₂sin(2ωt)+...+ A₆cos(6ωt)+ B₆sin(6ωt)

t=0:360;
x=a0+a'*cos(w*k'*t)+b'*sin(w*k'*t);
plot(t,x,'b')
hold off

Ahora t es un vector fila de 361 elementos, pero el producto w*k'*t es una matriz de 6×361 elementos. a o b son vectores columna de 6 elementos, dimensión 6×1. El producto a'*cos(w*k'*t) tiene la siguiente dimensión (1×6)×(6×361)=1×361, es decir, un vector fila de 361 elementos.

Unimos las porciones de código en un único script

t=0:30:330; 
x=[408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
hold on
plot(t,x,'ro','markersize',4,'markerfacecolor','r')

n=length(x);
w=2*pi/360;
a0=sum(x)/n;
k=1:6;
a=cos(w*k'*t)*x'*2/n;
b=sin(w*k'*t)*x'*2/n;

t=0:360;
x=a0+a'*cos(w*k'*t)+b'*sin(w*k'*t);
plot(t,x,'b')

xlim([0 360])
xlabel('Ascensión (grados)')
ylabel('Declinación (minutos)')
title('Posición del asteroide Pallas')
grid on
hold offf

Ejemplos del curso de Física

Funciones básicas: sqrt, sin, etc. y programación

Dinámica de una partícula unida a una goma elástica

Secuencia de colisiones elásticas (II)

Trayectorias hiperbólicas

Trayectoria de un proyectil disparado desde la superficie de la Tierra

Trayectoria de un proyectil disparado desde una altura h sobre la superficie de la Tierra

Desviación hacia el este de un cuerpo que cae

Choque de un meteorito con la Tierra

Encuentro de una sonda espacial con el planeta Júpiter

ondas/paraxial/paraxial.html

Macroestado, microestados. Temperatura, entropía

Vaciado de un depósito cerrado

Esfera cargada próxima a un plano conductor a potencial cero

Dos esferas conductoras una de las cuales está a potencial cero

Fuerza entre dos esferas conductoras

Ping-Pong eléctrico

Circuitos de corriente alterna

Intensidad de la radiación solar sobre una superficie plana inclinada

Producto escalar: dot, norm

Movimiento curvilíneo

Suma de los términos de una serie: sum

Sucesivos rebotes en un plano inclinado

Formulación discreta de las ecuaciones del movimiento de un cohete

Análisis de las alturas y periodos de las olas

Máximo o mínimo valor de un conjunto de datos: max, min

Respuesta de un oscilador a una fuerza impulsiva (I)

El arco iris

Ondas térmicas

Cohete propulsado por agua