Triangulación y convexidad

Este es un programa en MATLAB, clickea en el axes y luego de varias veces haz clic derecho. Verás la triangulación por Delaunay y la capsula convexa. Si desean ver un poco de teoría sobre convexidad y cápsula convexa, revisen : Functional Analisys with Applications , Zeidler.

clear all
clc
figure(1)
xlim([0 1])
ylim([0 1])
hold on
i=0;
X=[];
Y=[];
while 1
    [x,y,buttom]=ginput(1);
    if buttom==3
        break
    end
    plot(x,y,'ok')
    i=i+1;
    X(i)=x;
    Y(i)=y;  
end

% si usas Octave solo cambia esto por lo último

%T=delaunay(X',Y');

%k = convhull(X',Y');

%triplot(T,X,Y)

%plot(X'(k),Y'(k), '-rs')

dt=DelaunayTri(X',Y');
k = convexHull(dt);
triplot(dt)
plot(dt.X(k,1),dt.X(k,2), 'rs')

Listas en Python

La lista es un tipo de colección ordenada.  Sería equivalente a lo que en otros lenguajes se conoce por arrays o vectores. Las listas pueden contener cualquier tipo de dato: números, cadenas, booleanos, ... y también listas. Veremos a continuación algunos métodos útiles de estos objetos (todo en Python es un objeto).

L.append(object) : Añade un objeto al final de la lista.

L=[5,6,2,8,6]
x=5
L.append(x)
print L

[5, 6, 2, 8, 6, 5]

L=[5,6,2,8,6]
x=[6,6,3]
L.append(x)
print L

[5, 6, 2, 8, 6, [6, 6, 3]]

L.count(value) : Devuelve el número de veces que se encontró value en la lista.

L=[5,6,2,8,6]
L.count(6)

2

L.extend(iterable) : Añade los elementos del iterable a la lista.

L=[5,6,2,8,6]
x=[0,1,0,1,1]
L.extend(x)
print L

[5, 6, 2, 8, 6, 0, 1, 0, 1, 1]

L.index(value) : Devuelve la posición en la que se encontro la primera ocurrencia de value.

L=[5,6,2,8,6]
L.index(6)

1

L=[5,6,2,6,8,6]
L.index(6,2,4)

3

L.insert(index,object): Inserta el objeto en la posición index.

L=[5,6,2,6,8,6]
L.insert(2,[5,5])
print L

[5, 6, [5, 5], 2, 6, 8, 6]

L.pop([index]) : Devuelve el valor en la posición index y lo elimina de la lista. Si no seespecifica la posición, se utiliza el último elemento de la lista.

L=[5,6,2,6,8,6]
L.pop(2)

2

L.remove(value) : Elimina la priera ocurrencia de value en la lista.

L=[5,6,2,6,8,6]
L.remove(6)
print L
      [5, 2, 6, 8, 6]

L.reverse() : Invierte la lista.

L=[5,6,2,6,8,6]
L.reverse()
print L

[6, 8, 6, 2, 6, 5]

L.sort(cmp=None,key=None,reverse=false):  Ordena la lista. Si se especifica cmp este debe ser una función que tome como parametro dos valores x e y de la lista y devuelva -1 si x es menor que y, 0 si son iguales y 1 si x es mayor que y.

L=[5,6,2,6,8,6]
L.sort()
print L

[2, 5, 6, 6, 6, 8]

L=[5,6,2,6,8,6]
L.sort(reverse=True)
print L

[8, 6, 6, 6, 5, 2]

def f(x):
   return x*(x-1)
L=[0.1,0.5,0.25,0.6,0.15]
L.sort(key=f)
print L

[0.500000000000000, 0.600000000000000, 0.250000000000000, 0.150000000000000, 0.100000000000000]

Algunas funciones...

min([5,6,2,6,8,6])

2

max([5,6,2,6,8,6])

8

L=[5,6,2,8,6]
x=[5,6,2,8,6]
cmp(L,x)

0

L=list((4,5,5,5))   #las listas son mutables, las tuplas no.
print L

[4, 5, 5, 5]

L=[5,6,2,8,6]
len(L)

5

L1=[5,2,3]
L=sum(L1)
print L
10


map(function, sequence[,secuence, …]):

def f(n,m):
   return n*m+m
L1=[5,2,3]
L2=[8,1,4]
L=map(f,L1,L2)
print L

[48, 3, 16]

filter(function,sequence):
def f(n):
   return n**2%2==0
L1=[5,2,3]
L=filter(f,L1)
print L

[2]

reduce(function,sequence[,initial]):

def f(n,m):
   return n*m
L1=[5,2,3]
L=reduce(f,L1)
print L

30

Método de Euler en Python usando Sage

El clásico método de Euler(en Python):
La clase...

class ForwardEuler:
    def __init__(self, f, dt):
        self.f, self.dt = f, dt
    
    def set_initial_condition(self, u0, t0=0):
        self.u = []        
        self.t = []
        self.u.append(float(u0))
        self.t.append(float(t0))
        self.k = 0 
    
    def solve(self, T):
        tnew = 0
        while tnew <= T:
            unew = self.advance()
            self.u.append(unew)
            tnew = self.t[-1] + self.dt
            self.t.append(tnew)
            self.k += 1
        return numpy.array(self.u), numpy.array(self.t)
    
    def advance(self):
        u, dt, f, k, t = \
            self.u, self.dt, self.f, self.k, self.t[-1]
        unew = u[k] + dt*f(u[k], t)
        return unew

Utilizando la clase:

def _f1(u, t):
    return 0.2*u*(1-u)
 
u0 = 0.1
dt = 0.1
T = 40
method = ForwardEuler(_f1, dt)
method.set_initial_condition(u0, 0)
u, t = method.solve(T)
  
g1=line([[t[i],u[i]] for i in range(len(u))],rgbcolor=(1/4,1/8,3/4))
g2 = text("Ec. Log.", (20, 1))
g = g1+g2
show(g, xmin=0, xmax=41, ymin=0, ymax=1)

Curva Dragón

Curva dragón, así se llama este fractal, como siempre, con los fractales la idea es simple: "hacer lo mismo varias veces". Este es el primero que me salío. La figura inicial son dos lados consecutivos de un cuadrado. Lo rotamos 90º, y a toda la figura resultante, lo volvemos a rotar. Así sucesivamente.

Esta es una variante.

Código:

def rot_point(x,p):    
    size=len(x)
    P=[]
    for i in range(size):        
        P.append([-x[i][1]+p[1]+p[0],x[i][0]-p[0]+p[1]])
    return P
def curva_dragon(N):
    DRA=[[-1,0],[0,0],[0,-1]]
    for i in range(1,N):
        tem=[]
        tem=rot_point(DRA[0:2^i],DRA[2^i])
        tem.reverse()
        DRA=DRA+tem       
    return DRA
D=curva_dragon(10)
line(D)

Series de Taylor como ecuación en diferencias

Sabemos que :

para todo x.

Computacionalmente podemos calcular la expresión anterior como una ecuación en diferencias.

Es muy práctico ver el problema a calcular como un ecuación en diferencias.

Código en Python:

def exp_diffeq(x, N):
    n = 1
    an_prev = 1.0 # a_0
    en_prev = 0.0 # e_0
    while n <= N:
        en = en_prev + an_prev
        an = x/n*an_prev
        en_prev = en
        an_prev = an
        n += 1
    return en

Exp=exp_diffeq(1,10)
Print Exp

2.71828152557319

Día Juliano en Fortran

Desarrollado para el curso de Geofísica. Calcula el día juliano para: TANNO, TMES, TDIA, THORA, TMINUTO.

program GREJD
integer TMES, TANNO, TDIA, JDIA
real THORA,TMINUTO
TMINUTO=25
THORA=5                
TDIA=28                
TMES=7            
TANNO=1821
THORA=(THORA+(TMINUTO/60))/24+0.5
if (TMES > 2) then
   TMES=TMES-3                 
else
   TMES=TMES+9
TANNO=TANNO-1
endif
JDIA = (TANNO / 4000) * 1460969;
TANNO = MOD(TANNO,4000);
JDIA=JDIA +(((TANNO / 100)*146097)/4) +((MOD(TANNO,100) * 1461) / 4)+(((153 * TMES) + 2) / 5) +TDIA +1721119;
write (*,*) 'DIA JULIANO:'
PRINT *,'PARTE ENTERA',JDIA
PRINT *,'PARTE DECIMAL',dble(THORA)
end

Gráficas con Sage

Ejemplos

L = [[cos(pi*i/100),sin(pi*i/100)] for i in range(200)]

p = polygon(L, rgbcolor=(1,1,0))

show(p)

L=[[-1+cos(pi*i/100)*(1+cos(pi*i/100)),2*sin(pi*i/100)*(1-cos(pi*i/100))] for i in range(200)]
p = polygon(L, rgbcolor=(1/8,3/4,1/2))
show(p)

L = [[cos(pi*i/100)^3,sin(pi*i/100)] for i in range(200)]
p = line(L, rgbcolor=(1/4,1/8,3/4))
t = text("a bulb", (-1.7, 0.5))
x = text("x axis", (2,-0.2))
y = text("y axis", (0.6,1.3))
g = p+t+x+y
show(g, xmin=-1.5, xmax=2, ymin=-1, ymax=1.3)

L = [[sin(5*pi*i/100)^2*cos(pi*i/100)^3,sin(5*pi*i/100)^2*sin(pi*i/100)] for i in range(200)]

p = polygon(L, rgbcolor=(1/3,1/2,3/5))

show(p)

v=[(i,floor(i)) for i in range(-5,5)]
plot_step_function(v, vertical_lines=False)

g(x)=1/(x-1)

plot(g,(x,-2,2),ymin=-10, ymax=10,detect_poles=True,color='red')+ line([(1,-10), (1,10)],color='green',linestyle='--')

g1(x)=tan(x)

plot(g1,(x,-1.5*pi,1.5*pi),ymin=-10, ymax=10,detect_poles=True,color='red')+ line([(1.5*pi,-10), (1.5*pi,10)],color='green',linestyle='--')+ line([(-1.5*pi,-10), (-1.5*pi,10)],color='green',linestyle='--')+ line([(0.5*pi,-10), (0.5*pi,10)],color='green',linestyle='--')+ line([(-0.5*pi,-10), (-0.5*pi,10)],color='green',linestyle='--')