Teselación e interpolación de Datos Dispersos

esta demostración describe envolturas convexas, teselaciones Delaunay, y diagramas Voronoi de 3 dimensiones. Tambien muestra cómo interpolar datos dispersos tridimencionales.

Contenidos

Envoltura Convexa

Muchas aplicaciones en ciencia, ingeniería, estadística, y matemáticas usan estructuras como envolturas convexas, teselaciones Delaunay y diagramas Voronoi para el análisis de los datos. MATLAB® le permite analizar geométricamente conjuntos de datos en cualquier dimensión.

Aquí, en 3 dimensiones, se muestra un conjunto de 50 puntos con con su envoltura convexa.

% Crear los Datos.
n = 50;
X = randn(n,3);
% Graficar los Puntos.
plot3(X(:,1),X(:,2),X(:,3),'ko','markerfacecolor','k');

% Calcular la envoltura convexa.
C = convhulln(X);
% Graficar la envoltura convexa.
hold on
for i = 1:size(C,1)
   j = C(i,[1 2 3 1]);
   patch(X(j,1),X(j,2),X(j,3),rand,'FaceAlpha',0.6);
end

% Modificar la perspectiva

view(3), axis equal off tight vis3d; camzoom(1.2)

colormap(spring)

rotate3d on

Cubo de Datos

Podemos crear un conjunto de datos X de los 8 vertices de un cubo incluido su centro.

X es una matriz 9-por-3 en donde cada fila es la coordenada tridimensional de un punto.

% Crear X.
X = zeros(8,3);
X([5:8,11,12,15,16,18,20,22,24]) = 1; % Esquinas.
X(9,:) = [0.5 0.5 0.5]; % Centro.
% visualizar X.
cla reset; hold on
d = [1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3];
plot3(X(d,1),X(d,2),X(d,3),'b:');
plot3(X(:,1),X(:,2),X(:,3),'b.','markersize',20);
t = text(X(:,1),X(:,2),X(:,3), num2str((1:9)'));
set(t,'VerticalAlignment','bottom','FontWeight','bold', 'FontSize',12);
view(3); axis equal tight off vis3d; camorbit(10,0);
rotate3d on

Calculando la Envoltura Convexa

La envoltura convexa de un conjunto de datos es la región convexa más pequeña que contiene el conjunto de datos. La envoltura convexa del conjunto de datos del cubo X puede se4r calculada por CONVHULLN.

Para este conjunto de datos X, la envoltura convexa tiene 12 facetas, cada una correspondiente a una fila en K y graficada encima con X. el cubo es transparente de modo que se pueden ver todas las facetas y los puntos de datos.

% Calcular la envoltura convexa.
tri = convhulln(X);
% Graficar los datos
cla reset;
plot3(X(:,1),X(:,2),X(:,3),'ko','markerfacecolor','k');

% Graficar la envoltura convexa.
for i = 1:size(tri,1)
  c = tri(i,[1 2 3 1]);
  patch(X(c,1),X(c,2),X(c,3),i,'FaceAlpha', 0.9);
end

% Modificar la vista.
view(3); axis equal tight off vis3d
rotate3d on

Teselación Delaunay

Una teselación Delaunay tridimensional es un conjunto de tetrahedros en donde ninguno de los puntos de datos estan contenidos en circunferencia alguna del tetraedro. La teselación Delaunay del conjunto de datos X puede ser calculada mediante DELAUNAYN.

Las 12 filas de T representan los 12 tetraedros que particionan el conjunto de datos X.

% Calcular la teselación Delaunay.
tri = delaunayn(X);
% Graficar los datos.
plot3(X(:,1),X(:,2),X(:,3),'ko','markerfacecolor','k');

% Graficar la teselación.
for i = 1:size(tri,1)
  y = tri(i,[1 1 1 2; 2 2 3 3; 3 4 4 4]);
  x1 = reshape(X(y,1),3,4);
  x2 = reshape(X(y,2),3,4);
  x3 = reshape(X(y,3),3,4);
  patch(x1,x2,x3,(1:4)*i,'FaceAlpha',0.8);
end

% Modificar la vista.
view(3); axis equal tight off vis3d; camorbit(10,0)
rotate3d on

Diagrama Voronoi

Un diagrama Voronoi particiona los datos espaciales en regiones poliédricas, con una region para cada punto de datos. Cualquier lugar dentro de una región esta más cercano a su punto de datos que cualquier otro en el conjunto. El diagrama Voronoi del conjunto de datos del cubo X puede ser calculado mediante VORONOIN.

V es el conjunto de vertices Voronoi. C representa el conjuunto de regiones Voronoi. Para nuestro conjunto de datos X, C tiene 9 regiones Voronoi. Aqui mostramos una región Voronoi, la región para el punto central del cubo.

% Calcular el diagrama Voronoi.
[c,v] = voronoin(X);
% Graficar los datos.
plot3(X(d,1),X(d,2),X(d,3),'b:.',X(9,1),X(9,2),X(9,3),'k.','markersize',20);

% Graficar el diagrama Voronoi.
nx = c(v{9},:);
tri = convhulln(nx);
for i = 1:size(tri,1)
   patch(nx(tri(i,:),1),nx(tri(i,:),2),nx(tri(i,:),3),rand,'FaceAlpha',0.8);
end

% Modificar la vista.
view(3); axis equal tight off vis3d; camzoom(1.5); camorbit(20,0)
rotate3d on

Dispersion Multidimensional

GRIDDATAN interpola datos dispersos multidimensionales. Usa DELAUNAYN para teselar los datos, y luego los interpola basado en la teselación. empezamos con un conjunto de datos de 500 puntos aleatorios en 3 dimensiones y calculamos los valores de una función start with a data set of 500 random points in 3 dimensions and compute the values of a function, el cuadrado de la distancia desde el origen, en cada uno de esos puntos.

% Crear los Datos.
n = 500;
X = 2*rand(n,3)-1;
v = sum(X.^2,2);
% Hacer un dibujo para mostrar cómo se define X.
cla reset; hold on
plot3([0.02, 0.47],[0.02,0.57],[0.02,0.57],'k-','linewidth',2);
plot3(0,0,0,'bo','markerfacecolor','b');
cube = zeros(8,3);
cube([5:8,11,12,15,16,18,20,22,24]) = 1; % Esquinas
cube(9,:) = [0.5 0.5 0.5]; % Centro.
plot3(cube(d,1),cube(d,2),cube(d,3),'r:');
plot3(0.5,0.6,0.6,'go','markerfacecolor','g');
text(0.02,-0.2,0,'El Origen','fontsize',12,'fontweight','bold');
text(0.55,0.6,0.6,'Un punto en X','fontsize',12,'fontweight','bold');
text(0.28,0.3,0.35,'sqrt(v)','fontsize',12,'fontweight','bold');
view(3); axis equal tight off vis3d;
camorbit(20,-10);
rotate3d on

Con GRIDDATAN podemos interpolar X y los valores v sobre una cuadricula X0 para obtener los valores de la función v0 sobre esta cuadrícula.

Los puntos negros son X y los puntos rojos son la cuadrícula X0.

% Cuadricular los datos.
d  = -0.8:0.2:0.8;
[x0,y0,z0] = meshgrid(d,d,d);
X0 = [x0(:) y0(:) z0(:)];
v0 = reshape(griddatan(X,v,X0),size(x0));
% Graficar los resultados.
cla reset; hold on;
plot3(X(:,1),X(:,2),X(:,3),'k+','markerfacecolor','k');
plot3(X0(:,1),X0(:,2),X0(:,3),'r.','markerfacecolor','r');
view(3); axis equal tight off vis3d; camzoom(1.6);
rotate3d on

Para visualizar la superficie en todos los puntos en donde la función toma un valor constante, podemos usar ISOSURFACE y ISONORMALS. Dado que la función es el cuadrado de la distancia desde el origen squared distance from the origin, la supèrficie en un valor constante es una esfera.

c = 0.6;  % valor constante
cla reset; hold on;
h = patch(isosurface(x0,y0,z0,v0,c),'FaceColor','red','EdgeColor','none');
isonormals(x0,y0,z0,v0,h);

view(3); axis equal tight off vis3d; camzoom(1.6); camlight; lighting phong
rotate3d on

Con más puntos de datos en X y una cuadrícula más densa X0, la esfera es mas definida pero el cálculo toma más tiempo.

Aqui hay una esfera precalculada generada usando 5000 puntos de datos en X y una distancia ebtre los puntos de cuadrícula de 0.05.

% Load saved results.
load qhulldemo
cla reset; hold on
d = -0.8:0.05:0.8;
[x0,y0,z0] = meshgrid(d,d,d);
h = patch(isosurface(x0,y0,z0,v0,0.6));
isonormals(x0,y0,z0,v0,h);
set(h,'FaceColor','red','EdgeColor','none');
view(3); axis equal tight off vis3d; camzoom(1.6)
camlight; lighting phong
rotate3d on

Ver archivo .m [Teselación e Interpolación de Datos Dispersos]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

Traducción アタラクシア [Ataraxiainc]

Diferencia Finita Laplaciana

El ejemplo muestra el cálculo y la representación de la diferencia finita laplaciana sobre un dominio en forma de L.

Contenidos

El Dominio

Para este Ejemplo, numerosos puntos NUMGRID dentro de un dominio en forma de L. La Función SPY es una herramienta muy usada para visualizar el patrón de elementos diferentes de cero en una matriz dada

R = 'L'; % Otras posibles formas incluyen S,N,C,D,A,H,B
% Generar y visualizar la malla.
n = 32;
G = numgrid(R,n);
spy(G)
title('Una Red de Diferencias Finitas')

% Mostrar una versión reducida como muestra.

g = numgrid(R,12)

g =

     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     6    11    16    21    26    36    46    56    66     0
     0     2     7    12    17    22    27    37    47    57    67     0
     0     3     8    13    18    23    28    38    48    58    68     0
     0     4     9    14    19    24    29    39    49    59    69     0
     0     5    10    15    20    25    30    40    50    60    70     0
     0     0     0     0     0     0    31    41    51    61    71     0
     0     0     0     0     0     0    32    42    52    62    72     0
     0     0     0     0     0     0    33    43    53    63    73     0
     0     0     0     0     0     0    34    44    54    64    74     0
     0     0     0     0     0     0    35    45    55    65    75     0
     0     0     0     0     0     0     0     0     0     0     0     0

La Discreta Laplaciana

Se usa DELSQ para generar la laplaciana discreta. La Función SPY da un sentido gráfico de la población de la matrix.

D = delsq(G);
spy(D)
title('La Laplaciana de 5 Puntos')

% Número de puntos interiores

N = sum(G(:)>0)

N =

   675

El Problema del Valor límite de Dirichlet

Finalmente, vamos a resolver el problema del valor límite de Dirichlet para el systema lineal esparcido. El problema es planteado como sigue:

  delsq(u) = 1 en el interior,
  u = 0 sobre el límite.

rhs = ones(N,1);
if (R == 'N') % Para la disección anidada, desactivar el ordenamiento de grado mínimo.
    spparms('autommd',0)
    u = D\rhs;
    spparms('autommd',1)
else
    u = D\rhs; % Esto se usa para R=='L' como en este ejemplo
end

La Solución

Trazar un mapa de la solución sobre la red y mostrarlo como un mapa de contorno.

U = G;
U(G>0) = full(u(G(G>0)));
clabel(contour(U));
prism
axis square ij

Ahora mostrar la solución como una gráfico de malla.

colormap((cool+1)/2);
mesh(U)
axis([0 n 0 n 0 max(max(U))])
axis square ij

Ver archivo .m [Diferencia Finita Laplaciana]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

Traducción アタラクシア [Ataraxiainc]

Demostración Gráfica de Eigenvalores y Valores Singulares EIGSHOW.

 
Contenidos
 
Demostración Gráfica de Eigenvalores y Valores Singulares EIGSHOW.  

 EIGSHOW presenta un experimento gráfico que expone el efecto sobre el
 círculo unitario del mapeo inducido por varias matrices 2 por 2.
 Un pushbutton (pulsador) permite la elección del modo "eig" o el modo
 "svd".

 En modo eig, el mouse puede ser usado para mover el vector x alrededor
 del círculo unitario.  La trayectoria resultante de A*x es trazada.  El
 objetivo es encontrar vectores x de modo que A*x es paralelo a x.  De
 modo que x es un eigenvector de A.  La longitud de A*x es el eigenvalor
 correspondiente.

 En modo svd, el mouse mueve dos vectores unitarios perpendiculares, x e
 y.
 Los A*x y A*y resultantes son trazados.  Cuando A*x es perpendicular a
 A*y entonces x e y son vectores singulares a derechas, A*x y A*y son
 múltiplos de vectores singulares izquierdos, y las longitudes de A*x y
 A*y son los correspondientes valores singulares.

 El título del gráfico es un menú de las matrices seleccionadas,
 incluyendo algunas con menos de dos eigenvectores reales. EIGSHOW(A)
 inserta A, la cual debe ser 2-por-2, en el menú.

 Hay algunas cuestiones por considerar:
    ¿Cuales matrices son singulares?
    ¿Cuales matrices tienen eigenvalores complejos?
    ¿Cuales matrices tienen doble eigenvalor?
    ¿Cúales matrices tienen eigenvalores iguales a valores singulares?
    ¿Cúales matrices tienen formas canónicas no diagonales de Jordan?

if nargin == 0;
   initialize
elseif arg == 0
   action
elseif arg < 0
   setmode(arg)
else
   initialize(arg);
end

%------------------

function initialize(arg)

if nargin == 0
   arg = 6;
end

if isequal(get(gcf,'tag'),'eigshow');
   h = get(gcf,'userdata');
   mats = h.mats;
else
   set(gcf,'numbertitle','off','menubar','none')
   h.svd = 0;
   mats = {
      '[5/4 0; 0 3/4]'
      '[5/4 0; 0 -3/4]'
      '[1 0; 0 1]'
      '[0 1; 1 0]'
      '[0 1; -1 0]'
      '[1 3; 4 2]/4'
      '[1 3; 2 4]/4'
      '[3 1; 4 2]/4'
      '[3 1; -2 4]/4'
      '[2 4; 2 4]/4'
      '[2 4; -1 -2]/4'
      '[6 4; -1 2]/4'
      'randn(2,2)'};
end

if all(size(arg)==1)
   if (arg < length(mats))
      mindex = arg;
      A = eval(mats{mindex});
   else
      A = randn(2,2);
      S = ['[' sprintf('%4.2f %4.2f; %4.2f %4.2f',A.') ']'];
      mindex = length(mats);
      mats = [mats(1:mindex-1); {S}; mats(mindex)];
   end
else
   A = arg;
   if ischar(A)
      S = A;
      A = eval(A);
   else
      S = ['[' sprintf('%4.2f %4.2f; %4.2f %4.2f',A.') ']'];
   end
   if any(size(A) ~= 2)
      error('MATLAB:eigshow:InputSizeIncorrect','La Matriz debe ser 2 por 2')
   end
   mats = [{S};  mats];
   mindex = 1;
end

clf
if h.svd
    t = 'svd / (eig)';
else
    t = 'eig / (svd)';
end
uicontrol(...
   'style','pushbutton', ...
   'units','normalized', ...
   'position',[.86 .60 .12 .06], ...
   'string',t, ...
   'value',h.svd, ...
   'callback','eigshow(-1)');
uicontrol(...
   'style','pushbutton', ...
   'units','normalized', ...
   'position',[.86 .50 .12 .06], ...
   'string','help', ...
   'callback','helpwin eigshow')
uicontrol(...
   'style','pushbutton', ...
   'units','normalized', ...
   'position',[.86 .40 .12 .06], ...
   'string','close', ...
   'callback','close(gcf)')
uicontrol(...
   'style','popup', ...
   'units','normalized', ...
   'position',[.28 .92 .48 .08], ...
   'string',mats, ...
   'tag','mats', ...
   'fontname','courier', ...
   'fontweight','bold', ...
   'fontsize',14, ...
   'value',mindex, ...
   'callback','eigshow(get(gco,''value''))');

s = 1.1*max(1,norm(A));
axis([-s s -s s])
axis square
xcolor = [0 .6 0];
Axcolor = [0 0 .8];
h.A = A;
h.mats = mats;
h.x = initv([1 0]','x',xcolor);
h.Ax = initv(A(:,1),'Ax',Axcolor);
if h.svd
   h.y = initv([0 1]','y',xcolor);
   h.Ay = initv(A(:,2),'Ay',Axcolor);
   xlabel('Hacer A*x perpendicular a A*y','fontweight','bold')
   set(gcf,'name','svdshow')
else
   xlabel('Hcer A*x paralela a x','fontweight','bold')
   set(gcf,'name','eigshow')
end
set(gcf,'tag','eigshow', ...
   'userdata',h, ...
   'windowbuttondownfcn', ...
      'eigshow(0); set(gcf,''windowbuttonmotionfcn'',''eigshow(0)'')', ...
   'windowbuttonupfcn', ...
      'set(gcf,''windowbuttonmotionfcn'','''')')

%------------------

function h = initv(v,t,color)
h.mark = line(v(1),v(2),'marker','.','erase','none','color',color);
h.line = line([0 v(1)],[0 v(2)],'erase','xor','color',color);
h.text = text(v(1)/2,v(2)/2,t,'fontsize',12,'erase','xor','color',color);

%------------------

function action
h = get(gcf,'userdata');
pt = get(gca,'currentpoint');
x = pt(1,1:2)';
x = x/norm(x);
movev(h.x,x);
A = h.A;
movev(h.Ax,A*x);
if h.svd
   y = [-x(2); x(1)];
   movev(h.y,y);
   movev(h.Ay,A*y);
end

%------------------

function movev(h,v)
set(h.mark,'xdata',v(1),'ydata',v(2));
set(h.line,'xdata',[0 v(1)],'ydata',[0 v(2)]);
set(h.text,'pos',v/2);

%------------------

function setmode(arg)
h = get(gcf,'userdata');
h.svd = ~h.svd;
set(gcf,'userdata',h)
initialize(get(findobj(gcf,'tag','mats'),'value'))

EIGSHOW

Ver archivo .m [Demostración Gráfica de Eigenvalores y Valores Singulares EIGSHOW]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

Traducción アタラクシア [Ataraxiainc]

Matrices Exponenciales

Para información sobre el cálculo de matrices exponenciales, ver "Diecinueve formas dudosas de calcular el exponencial de una matriz, venticinco años despues" Revista SIAM 45, 3-49, 2003.

El portal Pseudospectra tambien es altamente recomendado. El sitio web es:

http://web.comlab.ox.ac.uk/projects/pseudospectra/

Aqui tenemos 3 de las 19 forma de calcular el exponencial de una matriz.

Contenidos

Iniciar con la Matriz A

A = [0 1 2; 0.5 0 1; 2 1 0]
Asave = A;

A =

         0    1.0000    2.0000
    0.5000         0    1.0000
    2.0000    1.0000         0

Escalando y Exponenciando

expmdemo1 es una implementación en código M del algoritmo 11.3.1 de Golub y Van Loan, Computación de Matrices, Tercera Edición.

% Escalar A por la potencia de 2 de modo que su norma sea < 1/2.
[f,e] = log2(norm(A,'inf'));
s = max(0,e+1);
A = A/2^s;

% aproximación Pade para exp(A)
X = A;
c = 1/2;
E = eye(size(A)) + c*A;
D = eye(size(A)) - c*A;
q = 6;
p = 1;
for k = 2:q
   c = c * (q-k+1) / (k*(2*q-k+1));
   X = A*X;
   cX = c*X;
   E = E + cX;
   if p
     D = D + cX;
   else
     D = D - cX;
   end
   p = ~p;
end
E = D\E;

% Deshacer el escalamiento mediante cuadratura repètida
for k = 1:s
    E = E*E;
end

E1 = E

Matriz Exponencial vía Series de Taylor

expmdemo2 usa la definición clásica para la matriz exponencial. como método numérico práctico, es lento a inadecuado si norm(A) es demasiado larga.

A = Asave;

% Series de Taylor para exp(A)
E = zeros(size(A));
F = eye(size(A));
k = 1;
while norm(E+F-E,1) > 0
   E = E + F;
   F = A*F/k;
   k = k+1;
end

E2 = E

E1 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

Matriz Exponencial vía Eigenvalues and Eigenvectors

expmdemo3 asume que la matriz tiene un conjunto completo de eigenvectores. Como un método numérico práctico practical numerical method, la precisión esta determinada por la condición de la matriz eigenvector.

A = Asave;

[V,D] = eig(A);
E = V * diag(exp(diag(D))) / V;

E3 = E

E3 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

Comparar los Resultados

Para esta matriz, todos lo hacen igualmente bien

E = expm(Asave);
err1 = E - E1
err2 = E - E2
err3 = E - E3

err1 =

  1.0e-014 *

    0.3553    0.1776         0
    0.0444    0.0444   -0.0444
   -0.0888   -0.0888   -0.3553
err2 =

  1.0e-014 *

         0         0   -0.2665
   -0.0888   -0.0888   -0.0888
    0.0888   -0.0888         0

err3 =

  1.0e-014 *

   -0.7994   -0.5329   -0.8882
   -0.7994   -0.6661   -0.8438
   -0.7994   -0.6217   -0.9770 

Fallos en las Series de Taylor

He aqui una matriz en donde los términos en las series de Taylor se hacen muy grandes antes de que lleguen a cero. Consecuentemente, expmdemo2 falla.

A = [-147 72; -192 93];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)

E1 =

   -0.0996    0.0747
   -0.1991    0.1494
E2 =

  1.0e+006 *

   -1.1985   -0.5908
   -2.7438   -2.0442

E3 =

   -0.0996    0.0747
   -0.1991    0.1494

Fallos de eigenvalres y Vectores

He aqui una matriz que no tiene un conjunto completo de eigenvectores. Consecuentemente, expmdemo3 fails.

A = [-1 1; 0 -1];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)

E1 =

    0.3679    0.3679
         0    0.3679
E2 =

    0.3679    0.3679
         0    0.3679

E3 =

    0.3679         0
         0    0.3679

 Ver archivo .m [Matrices Exponenciales]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

Traducción アタラクシア [Ataraxiainc]

Representación Gráfica de Matrices tipo Sparse

Este ejemplo muestra la malla de elementos finitos para un plano aerodinámico de la NASA, incluye dos alerones de arrastre.

Contenidos

Los datos son almacenados en el archivo AIRFOIL.MAT. Posee 4253 pares de coordenadas (x,y) de puntos de la malla. Tambien posee un arreglo de 12,289 pares de índices, (i,j), para especificar conexiones entre los puntos de malla.

load airfoil

La Malla de Elementos Finitos

Primero, escalar x e y por 2^(-32) para encajarlos en el rango [0,1]. Entonces formar la matriz adyacente de tipo sparse y hacerla definida positiva.

% Escalando x e y
x = pow2(x,-32);
y = pow2(y,-32);

% Formando la matriz adyacente de tipo sparse y haciéndola positiva
% definida.
n = max(max(i),max(j));
A = sparse(i,j,-1,n,n);
A = A + A';
d = abs(sum(A)) + 1;
A = A + diag(sparse(d));

% Graficando la malla de elementos finitos
gplot(A,[x y])

Visualizando el Patrón Sparse

SPY es usado para visualizar el patron sparse. SPY(A) grafica el patrón sparse de la matriz A.

spy(A)
title('La Matriz Adyacente.')

Reorganización Simétrica - Inversa Cuthill-McKee

SYMRCM usa la tecnica inversa Cuthill-McKee para reorganizar la matriz adyacente. r = SYMRCM(A) da como resultado un vector permutación r tal que A(r,r) tiende a tener sus elementos diagonales más cercanos a la diagonal que A. Este es un buen preordenamiento para la factorización LU o Cholesky que vienen de problemas "largos, cenceños" (...). funciona tanto para A simétrica como A asimétrica.

r = symrcm(A);
spy(A(r,r))
title('Inversa Cuthill-McKee')

Reorganización Simétrica - COLPERM

Use j = COLPERM(A) para obtener un vector permutación que reorganiza la columna de la matriz sparse A en orden no-decrciente de la cuenta de elementos diferentes de cero. Es algunas veces útil como preordenamiento de la factorización LU: lu(A(:,j)).

j = colperm(A);
spy(A(j,j))
title('Reorganización Cuenta Columna')

Reorganización simétrica - SYMAMD

SYMAMD da una permutación simétrica aproximada de minimo grado. p = SYMAMD(S), para una matriz definida simétrica positiva A, retorna el vector permutación p tal que S(p,p) tiende a tener un factor Cholesky más esparcido que S. Algunas veces SYMAMD funciona bien para matrices simétricas indefinidas tambien.

m = symamd(A);
spy(A(m,m))
title('Aproximado de Mínimo grado')

 Ver archivo ,m [Representación Gráfica de Matrices de Tipo Sparse]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

Traducción アタラクシア [Ataraxiainc]

Matrices de Tipo SPARSE

Este ejemplo muestra que reorganizar las filas y columnas de una matriz sparse S puede afectar el tiempo y almacenamiento requerido para una operación de matrices como la factorización de S en su descomposición Cholesky, S=L*L'

Contenidos

Visualizando una Matriz de Tipo SPARSE

Un gráfico SPY muestra los elementos diferentes de cero en una matriz.

Este gráfico spy muestra una matriz positiva simétrica definida SPARSE derivada de una porción de la matriz de prueba Harwell-Boeing "west0479", una matriz que describe las conecciones en un modelo de columna de difracción en una planta química.

load('west0479.mat')
A = west0479;
S = A * A' + speye(size(A));
pct = 100 / prod(size(A));

clf; spy(S), title('Una Matriz Simétrica de tipo Sparse')
nz = nnz(S);
xlabel(sprintf('Elementos diferentes de Cero=%d (%.3f%%)',nz,nz*pct));

Calculando el Factor Cholesky

Ahora vamos a calcular el factor Cholesky L, en donde S=L*L'. Notese que L contiene muchos más elementos diferentes a cero que S no factorizado, pues el cálculo de la factorización Cholesky crea un "fill-in" (Llenado) de elementos diferentes de cero. Esto ralentiza el algoritmo e incrementa el costo de almacenamiento.

tic, L = chol(S)'; t(1) = toc;
spy(L), title('Descomposición Cholesky de S')
nc(1) = nnz(L);
xlabel(sprintf('Elementos diferentes de cero=%d (%.2f%%)   Duración=%.2f sec',nc(1),nc(1)*pct,t(1)));

Reorganización para Acelerar el Cálculo

Al reordenar las filas y columnas de una matriz, sería posible reducir la cantidad de fill-in creado por la factorización, reduciendo así la duración y el costo de almacenamiento.

Vamos a probar tres diferentes organizaciones soportadas por MATLAB®.

Usando la Inversión Cuthill-McKee

El comando SYMRCM usa la Inversión Cuthill-McKee reorganizando el algoritmo para mover todos los elementos diferentes de cero cerca a la diagonal, reduciendo el "ancho de banda" de la matriz original.

p = symrcm(S);
spy(S(p,p)), title('S(p,p) Despues del ordenamiento Cuthill-McKee')
nz = nnz(S);
xlabel(sprintf('Elementos diferentes de Cero=%d (%.3f%%)',nz,nz*pct));

El fill-in producido por lan factorización Cholesky es confinado a la banda, de modo que esa factorización de la matriz reordenada toma menor duración y almacenamiento.

tic, L = chol(S(p,p))'; t(2) = toc;
spy(L), title('chol(S(p,p)) Tras el ordenamiento Cuthill-McKee')
nc(2) = nnz(L);
xlabel(sprintf('Elementos diferentes de Cero=%d (%.2f%%)   Duración=%.2f sec', nc(2),nc(2)*pct,t(2)));

Usando Cuenta Columna

El comando COLPERM usa el algoritmo de reorganización de columnas para mover las filas y columnas con superiores cuentas diferentes a cero hacia el final de la matriz.

q = colperm(S);
spy(S(q,q)), title('S(q,q) Tras el ordenamiento Cuenta Columna')
nz = nnz(S);
xlabel(sprintf('Elementos diferentes de Cero=%d (%.3f%%)',nz,nz*pct));

Para este ejemplo, el ordenamiento cuenta columna se hace para reducir la duración y almacenamiento de la factorización Cholesky, pero este comportamiento no puede ser el esperado en general.

tic, L = chol(S(q,q))'; t(3) = toc;
spy(L), title('chol(S(q,q)) Tras el ordenamiento Cuenta Columna')
nc(3) = nnz(L);
xlabel(sprintf('Elementos diferentes de Cero=%d (%.2f%%)   Duración=%.2f sec',nc(3),nc(3)*pct,t(3)));

Usando el Grado Mínimo

El comando SYMAMD usa el algoritmo del mínimo grado aproximado (una poderosa técnica grafo-teorética) graph-theoretic technique) para producir bloques largos de ceros en la matriz.

r = symamd(S);
spy(S(r,r)), title('S(r,r) Tras la organización de Mínimo Grado')
nz = nnz(S);
xlabel(sprintf('Elementos diferentes de Cero=%d (%.3f%%)',nz,nz*pct));

Los bloques de ceros producidos por el algoritmo de mínimo grado son almacenados durante la factorización Cholesky. Esto puede reducir significativamente la duración y costos de almacenamiento.

tic, L = chol(S(r,r))'; t(4) = toc;
spy(L), title('chol(S(r,r)) Tras la organización de Mínimo Grado')
nc(4) = nnz(L);
xlabel(sprintf('Elementos diferentes de Cero=%d (%.2f%%)   Duración=%.2f sec',nc(4),nc(4)*pct,t(4)));

Resumiendo los Resultados

labels={'Original','Cuthill-McKee','Cuenta Columna','Grado Mínimo'};

subplot(2,1,1)
bar(nc*pct)
title('Elementos diferentes de Cero tras la factorización Cholesky')
ylabel('Porcentaje');
set(gca,'xticklabel',labels)

subplot(2,1,2)
bar(t)
title('Duración para completar la factorización Cholesky')
ylabel('Segundos');
set(gca,'xticklabel',labels)

 Ver archivo .m [Matrices de tipo SPARSE]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

Traducción アタラクシア [Ataraxiainc]

Grafos y Matrices

Este ejemplo da una explicación de la relación entre gráfos y matrices y una buena aplicacion de las matrices tipo SPARSE.

Un grafo es un conjunto de nodos con conexiones específicas entre si. Un ejemplo es el grafo de conectividad del domo geodésico Buckminster Fuller geodesic dome (una molécula de carbono-60 parecida a un balón de fútbol).

En MATLAB®, el grafo del domo geodésico puede ser generado mediante la función BUCKY.

% Definir las variabbles.
[B,V] = bucky;
H = sparse(60,60);
k = 31:60;
H(k,k) = B(k,k);

% Visualizar las variables.
gplot(B-H,V,'b-');
hold on
gplot(H,V,'r-');
hold off
axis off equal

Un grafo puede ser representado mediante su matriz adyacencia.

Para construir la matriz adyacencia, los nodos son numerados desde 1 hasta N. Entonces el elemento (i,j) de la matriz es ajustado a 1 si el nodo i es conectado al nodo j, y 0 de otra forma.

% Defina una matriz A.
A = [0 1 1 0 ; 1 0 0 1 ; 1 0 0 1 ; 0 1 1 0];

% Haga un dibujo que muestre los nodos conectados.
cla
subplot(1,2,1);
gplot(A,[0 1;1 1;0 0;1 0],'.-');
text([-0.2, 1.2 -0.2, 1.2],[1.2, 1.2, -.2, -.2],('1234')', ...
    'HorizontalAlignment','center')
axis([-1 2 -1 2],'off')

% Haga un dibujo que muestre la matriz adyacencia.
subplot(1,2,2);
xtemp=repmat(1:4,1,4);
ytemp=reshape(repmat(1:4,4,1),16,1)';
text(xtemp-.5,ytemp-.5,char('0'+A(:)),'HorizontalAlignment','center');
line([.25 0 0 .25 NaN 3.75 4 4 3.75],[0 0 4 4 NaN 0 0 4 4])
axis off tight

Aquí estan los nodos en un hemisferio de la bucky-bola, numerados polígono por polígono.

subplot(1,1,1);
gplot(B(1:30,1:30),V(1:30,:),'b-');
for j = 1:30,
    text(V(j,1),V(j,2),int2str(j),'FontSize',10);
end
axis off equal

Para visualizar la matriz adyacencia de este hemisferio, usamos la función SPY para graficar la silueta de los elementos diferentes a cero.

Note que la matriz es simétrica, mientras el nodo i esta conectado al nodo j, el nodo j esta conectado al nodo i.

spy(B(1:30,1:30))
title('spy(B(1:30,1:30))')

Ahora vamos a extender nuestro esquema numérico a todo el grafo reflejando la numeración de un hemisferio en el otro.

[B,V] = bucky;
H = sparse(60,60);
k = 31:60;
H(k,k) = B(k,k);
gplot(B-H,V,'b-');
hold on
gplot(H,V,'r-');
for j = 31:60
    text(V(j,1),V(j,2),int2str(j), ...
        'FontSize',10,'HorizontalAlignment','center');
end
hold off
axis off equal

Finalmente, aqui tenemos un gráfico SPY de la matriz final de tipo sparse.

spy(B)
title('spy(B)')

En muchos grafos útiles, cada nodo esta conectado tan solo a otros pocos nodos. Como resultado, las matrices adyacencia contienen solo algunas pocas entradas distintas a cero por fila.

Este ejemplo ha mostrado una área en donde las matrices SPARSE se vuelven útiles.

gplot(B-H,V,'b-');
axis off equal
hold on
gplot(H,V,'r-');
hold off

Ver archivo .m [Grafos y Matrices]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7
アタラクシア

Matemática de Precisión Simple

Se daran algunos ejemplos de manejo aritmético y álgebra lineal con datos de precisión simple. Tam bien se muestra un ejemplo de M-File en donde los resultados se calculan adecuadamente como de precision sencilla o doble dependiendo de los datos ingresados.

Contenidos

Crear Datos de Precisión Doble

Primero vamos a crear algunos datos, los cuales son de precisión doble por defecto.

Ad = [1 2 0; 2 5 -1; 4 10 -1]

Ad =

     1     2     0
     2     5    -1
     4    10    -1

Convertir a precisión Simple

Podemos convertir datos a precisión simple mediante la función single.

A = single(Ad); % ó A = cast(Ad,'single');

Creando Ceros y Unos de Precisión Simple

Tambien podemos crear ceros y unos de precisión simple con sus respectivas funciones.

n=1000;
Z=zeros(n,1,'single');
O=ones(n,1,'single');

Vamos a ver las variables en el espacio de trabajo.

whos A Ad O Z n

  Name         Size            Bytes  Class     Attributes

  A            3x3                36  single              
  Ad           3x3                72  double              
  O         1000x1              4000  single              
  Z         1000x1              4000  single              
  n            1x1                 8  double              

Podemos ver que algunas de las variables son de tipo single y que la variable A (la version de precisión simple de Ad) toma la mitad del número de bytes de memoria para ser almacenada puesto que los datos simples solo requieren 4 bytes (32-bits), en tanto que los dobles necesitan 8 bytes (64-bits).

Aritmética y Álgebra Lineal

Podemos manejar aritmética estandar y álgebra lineal sobre simples.

B = A'    % Matriz traspuesta

B =

     1     2     4
     2     5    10
     0    -1    -1

whos B

  Name      Size            Bytes  Class     Attributes

  B         3x3                36  single 

Vemos los resultados de esta operación, B, es un simple.

C = A * B % Multiplicación de la matriz

C =

     5    12    24
    12    30    59
    24    59   117

C = A .* B % Aritmética element-wise, de elemento racional

C =

     1     4     0
     4    25   -10
     0   -10     1

X = inv(A) % Matriz inversa

X =

     5     2    -2
    -2    -1     1
     0    -2     1

I = inv(A) * A % Confirmar que el resultado es la matriz identidad

I =

     1     0     0
     0     1     0
     0     0     1

I = A \ A  % Mejor forma para hacer una división de matrices que inv

I =

     1     0     0
     0     1     0
     0     0     1

E = eig(A) % Eigenvalores

E =

    3.7321
    0.2679
    1.0000

F = fft(A(:,1)) % FFT

F =

   7.0000          
  -2.0000 + 1.7321i
  -2.0000 - 1.7321i

S = svd(A) % Descomposición de un valor singular

S =

   12.3171
    0.5149
    0.1577

P = round(poly(A)) % El polinomio característico de una matriz

P =

     1    -5     5    -1

R = roots(P) % Raices de un polinomio

R =

    3.7321
    1.0000
    0.2679

Q = conv(P,P) % Convolver dos vectores
R = conv(P,Q)

Q =

     1   -10    35   -52    35   -10     1

R =

stem(R); % Graficar el resultado

Un programa que funciona para cualquier precisión doble o simple

Ahora vamos a observar una función pra calcular bastantes terminos en la secuencia de Fibonacci de modo que el radio es menor que la máquina correcta epsilon (eps) para datos de tipo doble o simple.

% Cuantos términos se necesitan para obtener resultados de precisión
% simple?
fibodemo('single')

% Cuantos términos se necesitan para obtener resultados de precisión doble?
fibodemo('double')

% Ahora vamos a ver el código involucrado.
type fibodemo

% Observe que hemos inicializado varias de nuestras variables, |fcurrent|,
% |fnext|, y |goldenMean|, con valores dependientes del tipo de dato
% entrante, y la tolerancia |tol| depende tambien del tipo de dato. la
% precisión simple requiere que calculemos menos términos que el cálculo
% equivalente de doble precisión.

ans =

    19

ans =

    41

function nterms = fibodemo(dtype)
%FIBODEMO Used by SINGLEMATH demo.
% Calculate number of terms in Fibonacci sequence.

% Copyright 1984-2004 The MathWorks, Inc. 
% $Revision: 1.1.4.1 $  $Date: 2004/03/22 23:53:55 $

fcurrent = ones(dtype);
fnext = fcurrent;
goldenMean = (ones(dtype)+sqrt(5))/2;
tol = eps(goldenMean);
nterms = 2;
while abs(fnext/fcurrent - goldenMean) >= tol
    nterms = nterms + 1;
    temp  = fnext;
    fnext = fnext + fcurrent;
    fcurrent = temp;
end

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

アタラクシア

Aritmética de Enteros

Se da un ejemplo del funcionamiento aritmético de datos enteros en señale e imágen.

Contenidos

Cargar Datos de Señal Entera

Cargar conjuntos de datos de medición que comprenden señales de cuatro instrumentos, usando 8 y 16-bit A-to-D's dando como resultado datos almacenados como int8, int16 y uint16. El tiempo es almacenado como uint16.

load integersignal

% Observe las variables
whos Signal1 Signal2 Signal3 Signal4 Time1

  Name            Size            Bytes  Class     Attributes

  Signal1      7550x1              7550  int8                
  Signal2      7550x1              7550  int8                
  Signal3      7550x1             15100  int16               
  Signal4      7550x1             15100  uint16              
  Time1        7550x1             15100  uint16              

Graficar los Datos

Primero vamos a graficar dos de las señales para ver los rangos de señal.

plot(Time1, Signal1, Time1, Signal2);
grid;
legend('Señal1','Señal2');

Podemos ver que los valores estan dentro del rango de -128 y 127, que es lo que podríamos esperar para int8. Es probable que estos valores necesiten ser escalados pàra calcular el valor físico actual que representa la señal e.g. Voltios

Procesado de Datos

Podemos utilizar aritmética estandar para enteros como +, -, *, y /. Supongamos que deseamos encontrar la suma de las señales 1 y 2.

SumSig = Signal1 + Signal2; % Aqui sumamos las señales enteras.

Podemos activar avisos que nos indiquen si la aritmética se satura. En el caso de int8, esto puede ocurrir si el resultado de una operación esta por fuera del rango -128 a 127. Vamos a activar los avisos con el comando intwarning y a re-ejecutar.

intwarning('on')
SumSig = Signal1 + Signal2; % Aquí substraemos la señales enteras de nuevo.

Warning: Out of range value or NaN computed in integer arithmetic.

Sin embargo, llevar a cabo este control ralentiza bastante la ejecución. Es mejor activarlo durante el desarrollo del algoritmo y apagarlo para la ejecución final. Ahora vamos a graficar la suma de las señales y ver en dónde se saturan.

intwarning('off')

cla;
plot(Time1, SumSig);
hold on;
Saturated = (SumSig == intmin('int8')) | (SumSig == intmax('int8')); % En cuentra en dónde se satura
plot(Time1(Saturated),SumSig(Saturated),'rd');grid;
hold off;

Los marcadores muestran en donde se ha saturado la señal.

Cargar Datos de Imagen Entera

A continuación vamos a observar la aritmética sobre algunos datos de imagen.

street1=imread('street1.jpg'); % Carga los Datos de imagen.
street2=imread('street2.jpg');
whos street1 street2

  Name           Size                Bytes  Class    Attributes

  street1      480x640x3            921600  uint8              
  street2      480x640x3            921600  uint8               

Podemos ver que las imágenes son a color de 24-bit, almacenada como tres planos de datos uint8.

Desplegar Imágenes

Display first image.

cla;
image(street1); % Despliega la imagen.
axis equal; axis off

Desplegar segunda imagen

image(street2); % Despliega imagen.
axis equal; axis off

Escalar una Imagen

Podemos escalar la imagen mediante una constante de doble presición pero mantener las imágenes almacenadas como enteros. Por ejemplo,

duller = 0.5 * street2; % Escala la imagen con una constante doble pero crea un entero
whos duller

  Name          Size                Bytes  Class    Attributes

  duller      480x640x3            921600  uint8              

subplot(1,2,1);
image(street2);
axis off equal tight
title('Original');  % Despliega la imagen

subplot(1,2,2);
image(duller);
axis off equal tight
title('Oscurecida');    % Despliega la imagen

Añadir las imágenes

Podemos añadir las dos imágenes de la calle juntas y graficar el resultado fantasmal.

combined = street1 + duller; % Añade imágenes |uint8|
subplot(1,1,1)
cla;
image(combined); % Despliega la imagen
title('Combinado');
axis equal; axis off

 Ver archivo .m [Aritmética de Enteros]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

アタラクシア

Ajuste Óptimo de Funciones No Lineales

Esta es una demostración del ajuste óptimo de una función no lineal a un conjunto de datos. Usa FMINSEARCH, una implementación del algoritmo Nelder-Mead simplex (búsqueda directa) para minimizar una función no lineal de varias variables.

Primero, crear algunos datos de ejemplo y graficarlos.

t = (0:.1:2)';
y = [5.8955 3.5639 2.5173 1.9790 1.8990 1.3938 1.1359 1.0096 1.0343 ...
     0.8435 0.6856 0.6100 0.5392 0.3946 0.3903 0.5474 0.3459 0.1370 ...
     0.2211 0.1704 0.2636]';
plot(t,y,'ro'); hold on; h = plot(t,y,'b'); hold off;
title('Datos de entrada'); ylim([0 6])

La meta es ajustar la siguiente función con dos parámetros lineales y dos parámetros no lineales a los datos:

   y =  C(1)*exp(-lambda(1)*t) + C(2)*exp(-lambda(2)*t)

Para ajustar esta función, se ha creado la función FITFUN. Dado el parámetro no lineal (lambda) y los datos (t e y), FITFUN calcula el error en el ajuste para esta ecuación y actuasliza la línea (h).

type fitfun

function err = fitfun(lambda,t,y)
%FITFUN Used by FITDEMO.
%   FITFUN(lambda,t,y) returns the error between the data and the values
%   computed by the current function of lambda.
%
%   FITFUN assumes a function of the form
%
%     y =  c(1)*exp(-lambda(1)*t) + ... + c(n)*exp(-lambda(n)*t)
%
%   with n linear parameters and n nonlinear parameters.

A = zeros(length(t),length(lambda));
for j = 1:length(lambda)
   A(:,j) = exp(-lambda(j)*t);
end
c = A\y;
z = A*c;
err = norm(z-y);

 Hacer una suposición para la estimación inicial de lambda (start) e invocar FMINSEARCH. Esto minimiza el error retornado por FITFUN ajustando lambda. Da como resultado el valor final de lambda. Usar una función de salida para graficar los ajustes intermedios.

start = [1;0];
% Usamos una función anónima para pasar parámetros adicionales t, y, h a la
% funcion de salida
outputFcn = @(x,optimvalues,state) fitoutputfun(x,optimvalues,state,t,y,h);
options = optimset('OutputFcn',outputFcn,'TolX',0.1);
estimated_lambda = fminsearch(@(x)fitfun(x,t,y),start,options)

estimated_lambda =

    3.5897
    0.0030

Ver archivo .m [Ajuste Óptimo de Funciones No lineales]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

アタラクシア

Predicción de la Población Estadounidense

Este ejemplo es mas antiguo que MATLAB®. Inició como un ejercicio en "Computer Methods for Mathematical Computations", de Forsythe, Malcolm y Moler, publicado por Prentice-Hall en 1977.

Ahora, MATLAB y Handle Graphics® hacen mucho más fácil variar los parámetros y ver los resultados, pero, los principios matemáticos subyacentes se mantienen inalterados. Se demuestra que el uso de polinómios de modesto grado incluso para predecir el futuro mediante extrapolación de datos es un negocio arriegado

$Revisión: 6.00.0.0 $ $Date: 2009/02/03 02:23:54 $

He aquí datos del censo de la población estadounidense desde 1900 a 2000.

% Time interval
t = (1900:10:2000)';
 
% Population
p = [75.995 91.972 105.711 123.203 131.669 ...
150.697 179.323 203.212 226.505 249.633 281.422]';
 
% Plot
plot(t,p,'bo');
axis([1900 2020 0 400]);
title(
'Population of the U.S. 1900-2000');
ylabel('Millions');

Cual es su predicción para la población en el año 2010?

p

p =

   75.9950
   91.9720
  105.7110
  123.2030
  131.6690
  150.6970
  179.3230
  203.2120
  226.5050
  249.6330
  281.4220

Vamos a ajustar los datos con un polinómio en t y a usarlo para extrapolarlo en t = 2010. Los coeficientes en el polinómio son obtenidos resolviendo un sistema de ecuaciones lineales que involucran una matriz Vandermonde de 11 por 11, cuyos elementos son potencias de escalas de tiempo, A(i,j) = s(i)^(n-j);

n = length(t);
s = (t-1950)/50;
A = zeros(n);
A(:,end) = 1;
for j = n-1:-1:1, A(:,j) = s .* A(:,j+1); end

Los coeficientes c para un polinómio de grado d que ajusta los datos p son obtenidos resolviendo un sistema lineal de ecuaciones que involucran la última columna d+1 de la matriz Vandermonde:

   A(:,n-d:n)*c ~= p

Si d es menor que 10, tenemos más ecuaciones que incógnitas y es apropiada una solución de mínimos cuadrados. Si d es igual a 10, la ecuación pude ser resuelta exactamente y el polinómio interpolaría los datos. En cualquier caso el sistema es resuelto con el operador barra de MATLAB. Aquí estan los coeficientes para el ajuste cúbico. 

c = A(:,n-3:n)\p

c =

    1.2629
   23.7261
  100.3659
  155.9043
 

Ahora evaluamos el polinómio para cada año desde 1900 hasta 2010 y graficamos los resultados.

v = (1900:2020)';
x = (v-1950)/50;
w = (2010-1950)/50;
y = polyval(c,x);
z = polyval(c,w);

hold on
plot(v,y,'k-');
plot(2010,z,'ks');
text(2010,z+15,num2str(z));
hold 
off

Comparar el ajuste cúbico con el cuártico. Nótese que los puntos extrapolados son muy diferentes

c = A(:,n-4:n)\p;
y = polyval(c,x);
z = polyval(c,w);

hold on
plot(v,y,
'k-');
plot(2010,z,'ks');
text(2010,z-15,num2str(z));
hold 
off

Como el grado aumente, la extrapolación se hace más irregular.

cla
plot(t,p,
'bo'); hold on; axis([1900 2020 0 400]);
colors = hsv(8); labels = {'Datos'};
for
 d = 1:8
[Q,R] = qr(A(:,n-d:n));
R = R(1:d+1,:); Q = Q(:,1:d+1);
c = R\(Q'*p); 
% Igual que c = A(:,n-d:n)\p;
y = polyval(c,x); 
z = polyval(c,11);
plot(v,y,
'color',colors(d,:));
labels{end+1} = ['Grado = ' int2str(d)];
end
legend(labels,2)

 Ver archivo. m [Predicción de la Población Estadounidense]

Copyright 1984-2009 The MathWorks, Inc.
Published with MATLAB® 7.7

Traducción por アタラクシア

• Configurando la Computación multihilos Mendiante el Menú de Preferencias
• Acerca de la Función de Ayuda Usada en este Ejemplo
• Midiendo las Mejoras en el Rendimiento para una Operación Simple
• Midiendo las Mejoras en el Rendimiento para Múltiples operaciones
• Graficando las mejoras en el Rendimiento para las Operaciones
• Resultados del Análisis de rendimiento y otros Temas
• Deshabilitando la Computación Multihilos