<jupyter_start><jupyter_text>Aca implementare EDP elipticas<jupyter_code># este es para eliptica pero con condiciones de contorno de dirichlet
# queda pendiente para las condiciones de newman
def elipticasEDP(ax0=0,ay0=0,nx=10,my=10,h=1,itmax=100,g=None,f=None,limcondi=None):
    '''
    Se resuelve la EDP :  laplaciano(u)  + f * u = g 
    Para eso llenamos la matriz u[i,j]  , para luego procedemos a refinarla usando metodo de gauss-seidel 
    Esta funcion recibe como parametros lo siguiente : 
    - ax0 , ay0 = son el borde superior izquierdo de la matriz . Seria como el u[0,0]
    - nx = seria el punto n-esimo de la division horizontal de la malla 
    - my = seria el punto m-esimo de la division vertical de la malla
    - h = seria el valor de paso , tanto para el eje "x" y "y" 
    - itmax = es el numero de iteracion para el metodo de gauss seidel
    - g,f = funciones de dos variable correspondientes al problema 
    - limcondi = son las condiciones de contorno que se conoce
    Lo que devuelve esta funcion la matriz "u" que resuelve el problema
    '''
    # reservamos espacio para una matriz de (my+1)*(nx+1)
    u = np.ndarray( shape=(my+1,nx+1), dtype = np.float64 )  
    ## llenamos los limites de los contornos . 
    for j in range(0,my+1):
        y = ay0 + j*h
        u[0,j] = limcondi(ax0,y) # parte izquierda 
        u[nx,j] = limcondi(ax0 + nx*h,y) # parte derecha 
    for i in range(0,nx+1):
        x = ax0 + i*h
        u[i,0] = limcondi(x,ay0) # parte inferior
        u[i,my] = limcondi(x,ay0 + my*h) # parte superior 
        
    for j in range(1,my):
        y = ay0 + j*h
        for i in range(1,nx):
            x = ax0 + i*h
            u[i,j] = 1 # con esta condicion iniciamos , luego refinamos    
    # refinamos la malla con gauss seidel  
    for k in range(0,itmax):
        for j in range(1,my):
            y = ay0 + j*h
            for i in range(1,nx):
                x = ax0 + i*h
                # aplicando metodo de gauss seidel 
                v = u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1]
                # g(x,y) no se pasa como argumentos ya que estan definidas como variable global 
                u[i,j] = (v-h**2*g(x,y)) / (4-h**2*f(x,y))
    
    return u
        
# implementemos un ejemplo para comprobar si nuestra funcion , funciona 
'''
El problema a resolver es laplaciano(u) + (-1/25)*u = 0 en un cuadrado unidad
donde las condiciones de contorno es u = q 
q = cosh(1/5x) + cosh(1/5*y)

Entonces tenemos 
f = -1/25
g = 0 
ax0 = ay0 = 0
h = 1/8 
nx =  8
my = 8 
itmax = 75
limcondi= cosh(1/5x) + cosh(1/5*y)
'''
def contorno(x,y):
    return np.float64( np.cosh(1/5*x) + np.cosh(1/5*y) )
    
    
def fprob(x,y):
    return np.float64(-1/25)

def gprob(x,y):
    return 0

matriz = elipticasEDP(ax0=0,ay0=0,nx=8,my=8,h=1/8,itmax=75,f=fprob,g=gprob,limcondi=contorno)
matriz # la matriz aparecera de la foorma que tu no la quieres , aparecera invertida<jupyter_output><empty_output><jupyter_text>Resolvemos el siguiente problema  Resolvemos el problema pero pensando que aislado es T=0<jupyter_code>'''
El problema a resolver es laplaciano(u) = 0 en un cuadrado unidad
donde las condiciones de contorno son las indicadas en la imagen

Entonces tenemos 
f = 0
g = 0 
ax0 = ay0 = 0
h = 0.25
nx =  4
my = 4
itmax = 75
'''
def contorno(x,y):
    if x==0 or y==0:
        resultado = 0
    else:
        if x==1:
            resultado = 100*y
        elif y == 1:
            resultado = 100*x
        else:
            resultado = 0             
    return resultado
         
def fprob(x,y):
    return 0

def gprob(x,y):
    return 0

matriz = elipticasEDP(ax0=0,ay0=0,nx=20,my=20,h=1/20,itmax=300,f=fprob,g=gprob,limcondi=contorno)
matriz # te aparecera de una forma diferente al que ves 
# invertimos 
matriz = matriz.T[::-1]# mapa de clalor
plt.figure(figsize=(6,6))
plt.imshow(matriz, cmap='hot', interpolation='nearest')
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Ahora si resolvemos el problema con las condiciones de newman<jupyter_code># este es para eliptica pero con condiciones de contorno de dirichlet derecha y arriba , newman -> abajo ,izquierda
# Esta funcion es elipticas con una ligera modificacion
def elipticasEDP_mod(ax0=0,ay0=0,nx=10,my=10,h=1,itmax=100,g=None,f=None,limcondi=None):
    '''
    Se resuelve la EDP :  laplaciano(u)  + f * u = g 
    Para eso llenamos la matriz u[i,j]  , para luego procedemos a refinarla usando metodo de gauss-seidel 
    Esta funcion recibe como parametros lo siguiente : 
    - ax0 , ay0 = son el borde superior izquierdo de la matriz . Seria como el u[0,0]
    - nx = seria el punto n-esimo de la division horizontal de la malla 
    - my = seria el punto m-esimo de la division vertical de la malla
    - h = seria el valor de paso , tanto para el eje "x" y "y" 
    - itmax = es el numero de iteracion para el metodo de gauss seidel
    - g,f = funciones de dos variable correspondientes al problema 
    - limcondi = son las condiciones de contorno que se conoce
    Lo que devuelve esta funcion la matriz "u" que resuelve el problema
    '''
    # reservamos espacio para una matriz de (my+1)*(nx+1)
    u = np.ndarray( shape=(my+1,nx+1), dtype = np.float64 )  
    ## llenamos los limites de los contornos . 
    for j in range(0,my+1):
        y = ay0 + j*h
        u[0,j] = 1 # parte izquierda /// esta parte sera reemplazada con condiciones de newman
        u[nx,j] = limcondi(ax0 + nx*h,y) # parte derecha 
    for i in range(0,nx):
        x = ax0 + i*h
        u[i,0] = 1 # parte inferior /// esta parte sera reemplazad con condiciones de newman
        u[i,my] = limcondi(x,ay0 + my*h) # parte superior 
        
    for j in range(1,my):
        y = ay0 + j*h
        for i in range(1,nx):
            x = ax0 + i*h
            u[i,j] = 1 # con esta condicion iniciamos , luego refinamos  
    
    # refinamos la malla con gauss seidel  
    for k in range(0,itmax):
        for j in range(0,my):
            y = ay0 + j*h
            for i in range(0,nx):
                x = ax0 + i*h
                if i == 0 or j==0 :
                    if i == 0 and j==0:
                        v = u[i+1,j] + u[i+1,j] + u[i,j+1] + u[i,j+1]
                    elif i == 0 and j != 0:
                        v = u[i+1,j] + u[i+1,j] + u[i,j+1] + u[i,j-1]
                    elif i != 0 and j == 0:
                        v = u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j+1]
                        
                else:
                    # aplicando metodo de gauss seidel 
                    v = u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1]
                # g(x,y) no se pasa como argumentos ya que estan definidas como variable global 
                u[i,j] = (v-h**2*g(x,y)) / (4-h**2*f(x,y))
    return u


# resolvemos el problema
'''
El problema a resolver es laplaciano(u) = 0 en un cuadrado unidad
donde las condiciones de contorno son las indicadas en la imagen

Entonces tenemos 
f = 0
g = 0 
ax0 = ay0 = 0
h = 0.25
nx =  4
my = 4
itmax = 75
'''
def contorno(x,y):
    if x==0 or y==0:
        resultado = 0
    else:
        if x==1:
            resultado = 100*y
        elif y == 1:
            resultado = 100*x
        else:
            resultado = 0             
    return resultado
         
def fprob(x,y):
    return 0

def gprob(x,y):
    return 0
# estudia mas !!!!
matriz = elipticasEDP_mod(ax0=0,ay0=0,nx=50,my=50,h=1/50,itmax=900,f=fprob,g=gprob,limcondi=contorno)
matriz # te aparecera de una forma diferente al que ves 
# invertimos 
matriz = matriz.T[::-1]# mapa de clalor
plt.figure(figsize=(7,7))
plt.imshow(matriz, cmap='hot', interpolation='nearest')
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Problema 2 de la PC<jupyter_code># haremos del problema de la PC 
'''
El problema a resolver es laplaciano(u) = -\rho/epsilon en un cuadrado unidad
donde las condiciones de contorno son las indicadas en la imagen

Entonces tenemos 
f = 0
g =  177,08 × 10−12Cos(3πx)Sin(2πy)
ax0 = ay0 = 0
h = 0.25
nx =  4
my = 4
itmax = 75
'''
lamda1 = 1
lamda2 = 1
def contorno(x,y):
    if x==0 :
        resultado = lamda1*y**2
    elif x==1:
        resultado = 1
    elif y==0:
        resultado = lamda2*x**3
    elif y==1:
        resultado = 1
    return resultado

         
def fprob(x,y):
    return 0

def gprob(x,y):
    epsilon = 8.854
    return np.float64(-177.08/epsilon*np.cos(3*np.pi*x)*np.sin(2*np.pi*y))

matriz = elipticasEDP(ax0=0,ay0=0,nx=4,my=4,h=1/4,itmax=100,f=fprob,g=gprob,limcondi=contorno)
matriz # te aparecera de una forma diferente al que ves 
matriz = matriz.T[::-1]# mapa de clalor
matriz
plt.imshow(matriz, cmap='hot', interpolation='nearest')
plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Aca implementare EDP parabolicas<jupyter_code># demostrare que puedo simular lo siguiente para hacer la simulacion de calor<jupyter_output><empty_output><jupyter_text>FTCS scheme [aqui](https://en.wikipedia.org/wiki/FTCS_scheme)<jupyter_code># funcion para el FTCS metodo 
def parabolicFTCS(k=0,h=0,f=None,x1=0,gx1=0,x2=1,gx2=0,t0=0,gt0=0,tfinal=100):
    '''
    Esta funcion resuelve problemas de una dimension de la forma diff2x_u + f(x)*diff1x_u = difft_u
    Para eso te pide los siguientes argumentos 
    - k = cantidad paso temporal 
    - h = cantidad de paso espacial
    - f = funcion de una variable "x" en la ecuacion diferencial
    - x1 = extremo izquierdo espacial
    - gx1 = condicion de extremo izquierdo u(x=x1,t) = gx1
    - x2 = extremo derecho espacial
    - gx2 = condicion de extremo derecho  u(x=x2,t) = gx2
    - t0 = tiempo inicial de la simulacion
    - gt0 = es la condicion inicial , es decir u(x,t=t0) = gt0
    - tfinal = hasta que tiempo final se esta haciendo la simulacion 
    
    Esta funcion retorna un array donde cada columna representa el estado u(x,t=tiempo) en un tiempo especifico
    Es decir las filas son como a "x" como columna es a "tiempo"
    
    '''
    # calculamos cuantas iteraciones vamos a hacer para llegar al tiempo final
    pasos_temporal = int((tfinal-t0)/k)
    # calculamos el segemnto entre x1 y x2 en cuanto sera dividio
    pasos_espacial = int((x2-x1)/h)
    
    ########## creamos la matriz donde haremos todos los calculos 
    # el +1 hace la cuenta de la condicion inicial (OJO)
    u = np.zeros( shape=(pasos_espacial+1,pasos_temporal+1), dtype = np.float64 ) 
    
    ####### LLenamos las condiciones de contorno
    s1 = k/h**2
    s2 = k/h
    # llenamos la condicion inicial temporal y espacial
    for i in range(1,pasos_espacial):
        x = x1 + i*h
        u[i,0] = gt0(x)

    for j in range(0,pasos_temporal+1):
        t = t0 + j*k
        u[0,j] = gx1(t)
        u[pasos_espacial,j] = gx2(t)
    
    ###### LLenamos todos los tiempos 
    for j in range(0,pasos_temporal):
        t = t0 + j*k
        for i in range(1,pasos_espacial):
            x = x1 + i*h
            u[i,j+1] = (s1+s2/2*f(x))*u[i+1,j] + (1-2*s1)*u[i,j] + (s1-s2/2*f(x))*u[i-1,j]
    return u<jupyter_output><empty_output><jupyter_text>$$\huge{ e^{\left(-\pi^2\cdot t\right)} sin(\pi\cdot x) }$$<jupyter_code># probemos nuestra funcion
# solucion analitica  u = np.e**(-np.pi**2*t) * np.sin(np.pi*x)
def fprob(x):
    return 0
def gx1prob(t):
    return 0
def gx2prob(t):
    return 0

def gt0prob(x):
    return np.float64(np.sin(np.pi*x))


matriz = parabolicFTCS(k=0.005,h=0.1,f=fprob,x1=0,gx1=gx1prob,x2=1,gx2=gx2prob,
                      t0=0,gt0=gt0prob,tfinal=0.01)
matriz
def analitico(x,t):
    return np.float64(np.e**(-np.pi**2*t) * np.sin(np.pi*x))
t0 = 0
x1 = 0
k = 0.005
h = 0.1
matrizprueba = np.ndarray( shape=(11,3), dtype = np.float64 ) 
for j in range(0,3):
    t = t0 + j*k
    for i in range(0,11):
        x = x1 + i*h
        matrizprueba[i,j] = analitico(x,t).round(6)
matrizprueba<jupyter_output><empty_output><jupyter_text>crank-nicholson scheme  full implicit scheme  Dufort-Frankel scheme    simulacion que verifique visualmente mi resultado  Tenemos que resolver el siguiente problema<jupyter_code>%matplotlib inline
from ipywidgets import interactive
import seaborn as sb
# funcion para el FTCS metodo 
# para resolver el problema se hizo un ligero cambio 
# para hallar la condicion de extrema derecha que usare 
def parabolicFTCS_mod(k=0,h=0,f=None,x1=0,gx1=0,x2=1,gx2=0,t0=0,gt0=0,tfinal=100):
    '''
    Esta funcion resuelve problemas de una dimension de la forma diff2x_u + f(x)*diff1x_u = difft_u
    Para eso te pide los siguientes argumentos 
    - k = cantidad paso temporal 
    - h = cantidad de paso espacial
    - f = funcion de una variable "x" en la ecuacion diferencial
    - x1 = extremo izquierdo espacial
    - gx1 = condicion de extremo izquierdo diff1t_u(x=x1,t) = gx1
    - x2 = extremo derecho espacial
    - gx2 = condicion de extremo derecho  u(x=x2,t) = gx2
    - t0 = tiempo inicial de la simulacion
    - gt0 = es la condicion inicial , es decir u(x,t=t0) = gt0
    . tfinal = hasta que tiempo final se esta haciendo la simulacion 
    
    //// es bueno aclarar que h , k ,x1 ,x2 ,t0, tfinal deben ser coherente , es decir para que no ocurra cosas raras
    //// en hacer los pasos temporales y pasos espaciales 
    
    Esta funcion retorna un array donde cada columna representa el estado u(x,t=tiempo) en un tiempo especifico
    Es decir las filas son como a "x" como columna es a "tiempo"
    
    '''
    # calculamos cuantas iteraciones vamos a hacer para llegar al tiempo final
    pasos_temporal = int((tfinal-t0)/k)
    # calculamos el segemnto entre x1 y x2 en cuanto sera dividio
    pasos_espacial = int((x2-x1)/h)
    
    ########## creamos la matriz donde haremos todos los calculos 
    # el +1 hace la cuenta de la condicion inicial (OJO)
    u = np.zeros( shape=(pasos_espacial+1,pasos_temporal+1), dtype = np.float64 ) 
    
    ####### LLenamos las condiciones de contorno
    s1 = k/h**2
    s2 = k/h
    u[0,0] = gx1(t0)
    # llenamos la condicion inicial temporal y espacial
    for i in range(1,pasos_espacial):
        x = x1 + i*h
        u[i,0] = gt0(x)
    for j in range(0,pasos_temporal+1):
        t = t0 + j*k
        # la modificacion se hizo aca --------- u[0,j] = gx1(t)
        u[pasos_espacial,j] = gx2(t)
    
    ###### LLenamos todos los tiempos 
    for j in range(0,pasos_temporal):
        t = t0 + j*k
        for i in range(1,pasos_espacial):
            x = x1 + i*h
            u[i,j+1] = (s1+s2/2*f(x))*u[i+1,j] + (1-2*s1)*u[i,j] + (s1-s2/2*f(x))*u[i-1,j]
        # se agrego lo siguiente , aca se uso un forward de segundo orden 
        u[0,j] = (2*k*gx1(t) + u[2,j] -4*u[1,j])/3 
    return u
        

# probemos nuestra funcion
# solucion analitica  u = np.e**(-np.pi**2*t) * np.sin(np.pi*x)
def fprob(x):
    return 1/x
def gx1prob(t):
    return 0
def gx2prob(t):
    return 1

def gt0prob(x):
    return 0

kprob = 0.005
hprob = 0.1
pasos_tempo = 200
pasos_espaciales_para_el_mapa_de_calor = 100
matriz = parabolicFTCS_mod(k=kprob,h=hprob,f=fprob,x1=0,gx1=gx1prob,x2=1,gx2=gx2prob,
                      t0=0,gt0=gt0prob,tfinal=kprob*pasos_tempo).T
# procedemos a graficar . 
matriz[100],matriz.shape , kprob/hprob**2
x = np.linspace(-1,1,pasos_espaciales_para_el_mapa_de_calor)
y = np.linspace(-1,1,pasos_espaciales_para_el_mapa_de_calor)

x,y = np.meshgrid(x,y)

def f(i,x,y,matrix):
    radio = np.sqrt(x**2+ y**2)
    numeros = np.round(radio/hprob)
    for k in range(len(numeros)):
        for j in range(len(numeros[k])):
            posicion = int(numeros[k,j])
            if posicion >= 1/hprob:
                radio[k,j] = matrix[i,int(1/hprob)]
            else:
                radio[k,j] = matrix[i,int(posicion)]
                
            
    return radio

def plotting(i=0):
    z = f(i,x,y,matriz)
    plt.figure(figsize=(9, 7))
    # recuerda que estas limitando el vmax , podrias cambiarlo 
    ax = sb.heatmap(z,xticklabels= 10 ,yticklabels= 10,cmap="gist_heat",vmin=0,vmax=1)
    plt.title(f"Mapa de calor {i*kprob} s")
    ax.invert_yaxis()
    

interactive_plot = interactive(plotting,i=(0,pasos_tempo,1))
interactive_plot<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12093773257256264, 'fraction_numerical': 0.036405358471843216, 'mean_word_length': 2.664772727272727, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '203729', 'license': 'MIT', 'n_tokens_mistral': 6739, 'n_tokens_neox': 6400, 'n_words': 1746, 'path': 'Curso_Metodos_Numericos_2020_I/codigos_despues_del_parcial/EDP.ipynb', 'repo_name': 'alonso121198/Algoritmos_aprendidos', 'size': '16124'}