martes, 22 de mayo de 2012

Octave


GNU Octave es un lenguaje de alto nivel destinado para el cálculo numérico. Como indica su nombre es parte de proyecto GNU. MATLAB es considerado su equivalente comercial. Provee una interfaz sencilla, orientada a la línea de comandos (consola), que permite la resolución de problemas numéricos, lineales y no lineales, además permite la ejecución de scripts y puede ser usado como lenguaje orientado al procesamiento por lotes.

Octave nació alrededor del año 1988, y fue concebido originalmente para ser usado en un curso de diseño de reactores químicos para los alumnos de Ingeniería Química de la Universidad de Texas y la Universidad de Wisconsin-Madison. Posteriormente en el año 1992, se decide extenderlo y comienza su desarrollo a cargo de John W. Eaton.1 La primera versión alpha fue lanzada el 4 de enero de 1993. Un año más tarde, el 17 de febrero de 1994 aparece la versión 1.0.

Octave posee una gran cantidad de herramientas que permiten resolver problemas de algebra lineal, cálculo de raíces de ecuaciones no lineales, integración de funciones ordinarias, manipulación de polinomios, integración de ecuaciones diferenciales ordinarias y ecuaciones diferenciales algebraicas. Sus funciones también se pueden extender mediante funciones definidas por el usuario escritas en el lenguaje propio de Octave o usando módulos dinámicamente cargados escritos en lenguajes como C, C++ y Fortran entre otros


Referencias

Autómatas celulares con Octave

Un Autómata celular es un sistema dinámico discreto adecuado para modelar sistemas reales que pueden describirse como una colección de objetos que interactúan localmente unos con otros. Se compone de una lattice cuyos elementos son conocidos como células (por lo que el espacio en un AC es discreto), un conjunto de estados para los valores que las células pueden tomar, conjuntos de vecindades entre cada célula de la lattice y células en su cercanía y una función local para asignarle un estado a cada célula conforme avanza el tiempo de manera discreta.

Célula. Elemento básico de un AC. La lattice que conforma al AC está compuesto por elementos conocidos como ‘células’.

Condición de frontera. Situación particular que se implementa en lattices finitas para considerar el comportamiento de las células cuyos elementos se encuentren en las fronteras de las mismas.

Configuración de un AC. Dadas las células de un AC en un tiempo t, la configuración es el conjunto de valores de dichas células en ese tiempo.

Conjunto de estados. Conjunto de elementos que dan valor a las células de una lattice.

Función global. Función que toma una configuración del AC en el tiempo t y, aplicando la función local, devuelve la configuración del AC en el tiempo t + 1.

Función local o Función de transición. Función que, por cada célula de un AC, toma los valores de los elementos en la vecindad de dicha célula para el tiempo t y devuelve el valor de esa célula en el tiempo t + 1. 




Referencias en Autómatas Celulares:


Autómata celular para un tejdo excitable en célula Purkinje con dominio toroidal

Conceptos:

Las células de Purkinje (o neuronas de Purkinje) son un tipo de neurona localizada en la corteza cerebelosa. Su nombre proviene de su descubridor, el de Purkinje, una capa de la corteza en el cerebelo situada entre la capa molecular y la capa granulosa. Las células de Purkinje envían proyecciones inhibidoras hacia el núcleo cerebelar profundo, y constituyen la única salida de toda la coordinación motriz en la corteza cerebelar. la celula de purkinje es una neurona de gran tamaño


En geometría el toroide es la superficie de revolución generada por una curva plana cerrada que gira alrededor de una recta exterior coplanaria (el eje de rotación situado en su mismo plano) con la que no se interseca. Su forma se corresponde con la superficie de los objetos que en el habla cotidiana se denominan donuts,argollas, anillos, aros o roscas. Cuando la curva cerrada es una circunferencia, la superficie se denomina toro.




% =========================================================================
% TEJIDO EXCITABLE EN DOMINIO TOROIDAL, CON CAPACIDAD DE PROPAGACION 
% ANISOTROPICA. 
% =========================================================================

clear all; close all;
% --------------------------------------------------DIMENSIONES DE LA TRAMA
N=40;  %-----------------------------------------------numero de renglones
M=40;  %------------------------------------------------numero de columnas 


iteraciones = 10000; %-------------------numero de ciclos para la simulacion

% ----------------------------------------------COORDENADAS DE ESTIMULACION
% en el centro
%      xESTIMULO=round(N/2);
%      yESTIMULO=round(M/2);
% al azar
        xESTIMULO=ceil(rand()*N); 
        yESTIMULO=ceil(rand()*M);


% -------------------------------------------------------------------------
%                       INICIALIZACIONES
% -------------------------------------------------------------------------

Enuevo = zeros(N,M)-90; %---------------todas las celulas en reposo (-90 V)
Eactual=zeros(N,M)-90;
RF = zeros(N,M);        % ----------arreglo de caracteristicas refractarias
CT = zeros(N,M);        % -------------arreglo con estimulaciones tempranas
CC = zeros(N,M);        % ---arreglo con influencia de celulas vecinas (VN)
g = 1;                %---------valor de referencia para la conductividad (Velocidad de propagación)

gSN = g*ones(N,M);       %--------------------------------conductividad S-N
gOE = g*ones(N,M);       %--------------------------------conductividad O-E
gNS = g*ones(N,M);       %--------------------------------conductividad N-S
gEO = g*ones(N,M);       %--------------------------------conductividad E-O

% =========================================================================
% INICIALIZACION DE LOS ESTADOS:
% =========================================================================

fase = 4*ones(N,M);     %---------todas las celulas inicializadas en reposo

% Para inicializacion desfasada: descomentar el siguiente bloque

%for i=1:N
% for j=1:M
%    num = floor(4*rand);
%    switch num
%      case 0, fase(i,j)=0; Eactual(i,j)=floor(rand*120-90);
%      case 1, fase(i,j)=1; Eactual(i,j)=floor(rand*20+10);
%      case 2, fase(i,j)=2; Eactual(i,j)=floor(rand*80-70);
%      case 3, fase(i,j)=4; Eactual(i,j)=floor(rand*20-90);
%    end
% end
%end     

% =========================================================================
% DESCOMENTAR PARA PROPAGACION ANISOTROPICA Y BLOQUEO
% -------------------------------------------------------------------------
%   xESTIMULO=round(N/4);
%   yESTIMULO=round(M/10);
%   gSN(1:round(N/2),round(M/3)) = zeros(round(N/2),1);
%   gOE(6:round(N/2)+5,round(M/3)) = zeros(round(N/2),1);
  
%   gNS(1:round(N/2),round(M/3)) = 5*g*ones(round(N/2),1);
%   gEO(1:round(N/2),round(M/3)) = 5*g*ones(round(N/2),1);
% -------------------------------------------------------------------------


% =========================================================================
% A PARTIR DE AQUI NO TOCAR!!!!!!!!!!!
% -------------------------------------------------------Excitacion inicial
estimulo = zeros(N,M);
estimulo(xESTIMULO,yESTIMULO) = 100;

figure 



% -------------------------------------------------------------------------
%                       CICLO PRINCIPAL
% -------------------------------------------------------------------------

for t = 1:iteraciones,

% -------------------------------------------------------------------------
%               CALCULO DE LAS INFLUENCIAS LOCALES (DOMINIO TOROIDAL)
% -------------------------------------------------------------------------
    for i=1:N,   
        for j=1:M,
           % -------------------------------------Esquinas de la rep. plana
           if (i==1) & (j==1),
              CC(i,j)=gSN(N,j)*(Eactual(N,j)-Eactual(i,j))+...
                               gOE(i,M)*(Eactual(i,M)-Eactual(i,j))+...
                               gNS(i+1,j)*(Eactual(i+1,j)-Eactual(i,j))+...
                               gEO(i,j+1)*(Eactual(i,j+1)-Eactual(i,j));
           elseif (i==1) & (j==M),
              CC(i,j)=gSN(N,j)*(Eactual(N,j)-Eactual(i,j))+...
                               gOE(i,j-1)*(Eactual(i,j-1)-Eactual(i,j))+...
                               gNS(i+1,j)*(Eactual(i+1,j)-Eactual(i,j))+...
                               gEO(i,1)*(Eactual(i,1)-Eactual(i,j));
           elseif (i==N) & (j==1),
              CC(i,j)=gSN(i-1,j)*(Eactual(i-1,j)-Eactual(i,j))+...
                               gOE(i,M)*(Eactual(i,M)-Eactual(i,j))+...
                               gNS(1,j)*(Eactual(1,j)-Eactual(i,j))+...
                               gEO(i,j+1)*(Eactual(i,j+1)-Eactual(i,j));
           elseif (i==N) & (j==M),
              CC(i,j)=gSN(i-1,j)*(Eactual(i-1,j)-Eactual(i,j))+...
                               gOE(i,j-1)*(Eactual(i,j-1)-Eactual(i,j))+...
                               gNS(1,j)*(Eactual(1,j)-Eactual(i,j))+...
                               gEO(i,1)*(Eactual(i,1)-Eactual(i,j));
           elseif (i==1)|(j==1)|(i==N)|(j==M),
               
           % ---------------------------------------Bordes de la rep. plana
              if i==1,
                 CC(i,j)=gSN(N,j)*(Eactual(N,j)-Eactual(i,j))+...
                               gOE(i,j-1)*(Eactual(i,j-1)-Eactual(i,j))+...
                               gNS(i+1,j)*(Eactual(i+1,j)-Eactual(i,j))+...
                               gEO(i,j+1)*(Eactual(i,j+1)-Eactual(i,j));
              end,   
              if j==1,
                 CC(i,j)=gSN(i-1,j)*(Eactual(i-1,j)-Eactual(i,j))+...
                               gOE(i,M)*(Eactual(i,M)-Eactual(i,j))+...
                               gNS(i+1,j)*(Eactual(i+1,j)-Eactual(i,j))+...
                               gEO(i,j+1)*(Eactual(i,j+1)-Eactual(i,j));
              end,   
              if i==N,
                 CC(i,j)=gSN(i-1,j)*(Eactual(i-1,j)-Eactual(i,j))+...
                               gOE(i,j-1)*(Eactual(i,j-1)-Eactual(i,j))+...
                               gNS(1,j)*(Eactual(1,j)-Eactual(i,j))+...
                               gEO(i,j+1)*(Eactual(i,j+1)-Eactual(i,j));
              end,   
              if j==M,
                 CC(i,j)=gSN(i-1,j)*(Eactual(i-1,j)-Eactual(i,j))+...
                               gOE(i,j-1)*(Eactual(i,j-1)-Eactual(i,j))+...
                               gNS(i+1,j)*(Eactual(i+1,j)-Eactual(i,j))+...
                               gEO(i,M)*(Eactual(i,M)-Eactual(i,j));
             end,   
           else
               
           % ----------------------------Puntos interiores de la rep. plana
             CC(i,j)=gSN(i-1,j)*(Eactual(i-1,j)-Eactual(i,j))+...
                              gOE(i,j-1)*(Eactual(i,j-1)-Eactual(i,j))+...
                              gNS(i+1,j)*(Eactual(i+1,j)-Eactual(i,j))+...
                              gEO(i,j+1)*(Eactual(i,j+1)-Eactual(i,j));
           end;
           
           % -----------------------------------Repeticion de la Excitacion
           if rem(t,5),
              CC(i,j) = CC(i,j) + estimulo(i,j);
           end;

        end;
    end;

% -------------------------------------------------------------------------
%               ACTUALIZACION DE POTENCIALES
% -------------------------------------------------------------------------
    for i=1:N,   
        for j=1:M,
           [fase(i,j),Enuevo(i,j),RF(i,j),CT(i,j)] = ...
         AC_celulaPURKINJE(fase(i,j),Eactual(i,j),RF(i,j),CT(i,j),CC(i,j));
           Eactual(i,j) = Enuevo(i,j);
        end;
    end;
    
    % ------------------------------------------Actualizacion de la grafica
    pcolor(Enuevo),colormap jet,shading flat;    
    axis square;
    drawnow;
    
end,
function [fase,Enuevo,RF,CT] = celulaPURKINJE(fase,Eactual,RF,CT,CC);
INH = 0.001;

switch fase,
    case 0,
        Enuevo = Eactual + 40;
        if Enuevo >= 30,
            fase = 1;
        end
    case 1,
        Enuevo = 0.7*Eactual;
        if Enuevo <= 10,
            fase = 2;
        end
    case 2,
        Enuevo = 1.1*Eactual + 4*CT - 6;
        if Enuevo <= -70,
            fase = 4;
            CT = 0;
            RF = 0;
        end
    case 4,
        Enuevo = Eactual + CC*RF;
        if Enuevo >= -60,
            fase = 0;
        else
            RF = min(RF + INH,1);
            CT = min(CT + INH,1);
        end
end

No hay comentarios:

Publicar un comentario