#!/usr/bin/python # -*- coding: iso-latin-1 -*- # Funciones para resolver Sistemas Dinámicos # # # # # ########################################### from math import * from random import * ##### SOLUCION DE SISTEMAS DE ECUACIONES DIFERENCIALES def Euler(X,F,dim,n,h=0.0001,arch='datos.dat',cada=1): """ Solución de un sistema de ecuaciones diferenciales de primer orden (al estilo de los Sistemas Dinámicos) Usando el método de Euler Parámetros: X, el vector inicial F, el de funciones del SD dim, la dimensión del SD n, número de puntos a evaluar h=0.0001, precisión (tamaño de paso) arch='datos.dat', el archivo de datos """ def write(X, t): for j in DIM: print>>ff,X[j], print>>ff, t X1=[0.0]*dim # arreglo auxiliar DIM = range(dim) # se guarda la condición inicial ff=open(arch,"w") for j in range(dim): print>>ff,X[j], print>>ff, 0 # método de Euler for t in range(n): for j in range(dim): X1[j] = X[j] + h*F[j](X) if t % cada == 0: write(X1, t+1) X, X1 = X1, X ff.close() def Euler2(X,F,dim,n,h=0.0001,arch='SDdatos.dat'): """ Notas: Solución de un sistema de ecuaciones diferenciales de primer orden (al estilo de los Sistemas Dinámicos) Usando el método de Euler de Segundo Orden: f'(x) = (f(x+h)-f(x-h)) / 2h Requiere dos puntos de inicio: Y_n+2 = Y_n + 2h F(X_n+1,Y_n+1) Parámetros: X, el vector inicial F, el de funciones del SD dim, la dimensión del SD n, número de puntos a evaluar h=0.0001, precisión (tamaño de paso) arch='datos.dat', el archivo de datos """ def write(X, t): for j in DIM: print>>ff,X[j], print>>ff, t X0=X[:] DIM = range(dim) # se guarda la condicin inicial ff=open(arch,"w") write(X0, 0) # evaluación del punto 2 X1 = [ X0[j] + h*h*F[j](X0) for j in DIM ] write(X1, 1) # método de Euler ORDEN 2 for t in range(n): X2 = [ X0[j] + 2*h*F[j](X1) for j in DIM ] write(X2, t+2) X0, X1 = X1, X2 del X2 ff.close() ### def RungeKutta(X,F,dim,n,h=0.0001,arch='SDdatos.dat',cada=10): """ Notas: Solución de un sistema de ecuaciones diferenciales de primer orden (al estilo de los Sistemas Dinámicos) Usando el método deRunge-Kutta de Cuarto Orden: Parámetros: X, el vector inicial F, el de funciones del SD dim, la dimensión del SD n, número de puntos a evaluar h=0.0001, precisión (tamaño de paso) arch='datos.dat', el archivo de datos """ def write(X, t): for j in DIM: print>>ff, X[j], print>>ff, t X0=X[:] DIM = range(dim) # se guarda la condicin inicial ff=open(arch,"w") write(X0, 0) # método de Runge-Kutta for t in range(n): k1 = [ h*F[j](X0) for j in DIM ] k1x = [ X0[i]+0.5*k1[i] for i in DIM ] k2 = [ h*F[j](k1x) for j in DIM ] k2x = [ X0[i]+0.5*k2[i] for i in DIM ] k3 = [ h*F[j](k2x) for j in DIM ] k3x = [ X0[i]+k3[i] for i in DIM ] k4 = [ h*F[j](k3x) for j in DIM ] X1 = [ X0[i] + 1./6.*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) for i in DIM ] if t % cada == 0: write(X1, t+1) del k1, k2, k3, k4, k1x, k2x, k3x X0, X1 = X1, X0 ff.close() ### def Poincare(X,F,L,v,dim,nL,h=0.0001,arch='SDPoincare'): """ Notas: Solución de un sistema de ecuaciones diferenciales de primer orden (al estilo de los Sistemas Dinámicos) Usando el método deRunge-Kutta de Cuarto Orden: Parámetros: X, el vector inicial F, el de funciones del SD L, lista de valores (planos) v, sobre qué variable dim, la dimensión del SD nL, número de puntos a evaluar h=0.0001, precisión (tamaño de paso) arch='datos.dat', el archivo de datos """ def write(X, t, s): for j in DIM: print>>f[s], X[j], print>>f[s], t X0=X[:] DIM = range(dim) # abrir archivos archant = arch f = [0]*len(L) for s in range(len(L)): arch = archant + "_" + str(L[s]) + ".dat" f[s] = open(arch,"w") # método de Runge-Kutta for t in range(nL): k1 = [ h*F[j](X0) for j in DIM ] k1x = [ X0[i]+0.5*k1[i] for i in DIM ] k2 = [ h*F[j](k1x) for j in DIM ] k2x = [ X0[i]+0.5*k2[i] for i in DIM ] k3 = [ h*F[j](k2x) for j in DIM ] k3x = [ X0[i]+k3[i] for i in DIM ] k4 = [ h*F[j](k3x) for j in DIM ] X1 = [ X0[i] + 1./6.*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) for i in DIM ] del k1, k2, k3, k4, k1x, k2x, k3x # check points at planes for s in range(len(L)): if abs(X0[v]-L[s]) < h: write(X1,t+1,s) X0, X1 = X1, X0 for s in range(len(L)): f[s].close() ### def ReaccionDifusionCA(X,F,A,D,dim,ter,n,CA,DimCA,h=0.0001, arch_datos='datos.dat', arch_mapa='map.dat'): """ Solución de un sistema de ecuaciones diferenciales de primer orden (al estilo de los Sistemas Dinámicos) con un término de difusión Usando el método de Euler de Segundo Orden: X' _0= A[0][0] X_1 + A[1] X_1 X_2 + ... Y_n+2 = Y_n + 2h F(X_n+1,Y_n+1) parámetros: X, punto de inicio F, conjunto de funciones de las variables A, matriz de coeficientes D, coeficientes de difusión dim, dimensión del SD ter, número máximo de términos en SD n, número de puntos a evaluar CA, matriz Autómata DimCA, tamaño de la matriz h, tamaño de paso (dt) arch, nombre del archivo de datos """ def write(X): for j in DIM: print>>ff,X[j], print>>ff X0 = X[:] DIM = range(dim) # se guarda la condicin inicial ff=open(arch,"w") write(X0) # evaluación del punto 2 usando Euler X1 = [ X0[j] + h**2*F[j](X0) for j in DIM ] write(X1) # método de Euler ORDEN 2 for t in range(n): X2 = [ X0[j] + 2*h*F[j](X1) for j in DIM ] write(X2) X0, X1 = X1, X2 ##### FUNCIONES APLICADAS A MATRICES (PENSANDO EN AUTOMATAS CELULARES) def Lap(x,i,j): Lx = x[i+1][j]+x[i-1][j]-2*x[i][j] Ly = x[i][j+1]+x[i][j-1]-2*x[i][j] return Lx+Ly ##### UN EJEMPLO ASQUEROSO ############################ X = [0.3, -0.23] A = [1.0, -1.0] D = [1.0, 1.0] DIM = 2 TER = 1 b = 0.0 DimCA = 100 CA = [] for k in range(DimCA): CA.append([random()]*DimCA) def F0(X): return A[0]*X[1] #+D[0]*Lap(CA[i][j]) def F1(X): return A[1]*X[0] #- b*X[1] #+D[1]*Lap(CA[i][j]) F = [F0, F1] ### primero probemos el caso SD: RungeKutta(X,F,DIM,50000,h=0.0001,arch='F0.dat') Euler2(X,F,DIM,50000,h=0.0001,arch='E0.dat') #ReaccionDifusionCA(X,F,A,D,2,1,10000,CA,DimCA,0.001,"RD.dat")