Categoría: Unidades

Unidades de estudio

  • 8.1 Regresión vs interpolación

    Referencia: Chapra 17.1 p 466. Burden 8.1 p498

    Cuando los datos de un experimento presentan variaciones o errores sustanciales respecto al modelo matemático, la interpolación polinomial presentada en la Unidad 4 es inapropiada para predecir valores intermedios.

    Minimos Cuadrados Lagrange 01

    En el ejemplo de Chapra 17.1 p470, se presentan los datos de un experimento mostrados en la imagen y la siguiente tabla:

    xi = [1,   2,   3,  4,  5,   6, 7]
    yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]

    Un polinomio de interpolación, por ejemplo de Lagrange de grado 6 pasará por todos los puntos, pero oscilando.

    Una función de aproximación que se ajuste a la tendencia general, que no necesariamente pasa por los puntos de muestra puede ser una mejor respuesta. Se busca una "curva" que minimice las diferencias entre los puntos y la curva, llamada regresión por mínimos cuadrados.


    Descripción con la media yi

    Considere una aproximación para la relación entre los puntos xi, yi como un polinomio, grado 0 que sería la media de yi. Para este caso, los errores se presentan en la gráfica:

    Minimos Cuadrados 03

    Otra forma sería aproximar el comportamiento de los datos es usar un polinomio de grado 1. En la gráfica se pueden observar que para los mismos puntos el error disminuye considerablemente.

    El polinomio de grado 1 recibe el nombre de regresión por mínimos cuadrados, que se desarrolla en la siguiente sección.

    Instrucciones Python

    # representación con la media
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    xi = [1,   2,   3,  4,  5,   6, 7]
    yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]
    
    # PROCEDIMIENTO
    xi = np.array(xi,dtype=float)
    yi = np.array(yi,dtype=float)
    n = len(xi)
    
    xm = np.mean(xi)
    ym = np.mean(yi)
    
    # SALIDA
    print('ymedia = ', ym)
    
    # grafica
    plt.plot(xi,yi,'o',label='(xi,yi)')
    plt.stem(xi,yi,bottom=ym, linefmt = '--')
    plt.xlabel('xi')
    plt.ylabel('yi')
    plt.legend()
    plt.show()
    

    Instrucciones Python - compara interpolación y regresión.

    Para ilustrar el asunto y para comparar los resultados se usa Python, tanto para interpolación y mínimos cuadrados usando las librerías disponibles. Luego se desarrolla el algoritmo paso a paso.

    # mínimos cuadrados
    import numpy as np
    import scipy.interpolate as sci
    import matplotlib.pyplot as plt
    
    # INGRESO
    xi = [1,   2,   3,  4,  5,   6, 7]
    yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]
    
    # PROCEDIMIENTO
    xi = np.array(xi)
    yi = np.array(yi)
    n = len(xi)
    
    # polinomio Lagrange
    px = sci.lagrange(xi,yi)
    xj = np.linspace(min(xi),max(xi),100)
    pj = px(xj)
    
    # mínimos cuadrados
    A = np.vstack([xi, np.ones(n)]).T
    [m0, b0] = np.linalg.lstsq(A, yi, rcond=None)[0]
    fx = lambda x: m0*(x)+b0
    fi = fx(xi)
    
    # ajusta límites
    ymin = np.min([np.min(pj),np.min(fi)])
    ymax = np.max([np.max(pj),np.max(fi)])
    
    # SALIDA
    plt.subplot(121)
    plt.plot(xi,yi,'o',label='(xi,yi)')
    plt.plot(xj,pj,label='P(x) Lagrange')
    plt.ylim(ymin,ymax)
    plt.xlabel('xi')
    plt.ylabel('yi')
    plt.title('Interpolación Lagrange')
    
    plt.subplot(122)
    plt.plot(xi,yi,'o',label='(xi,yi)')
    plt.plot(xi,fi,label='f(x)')
    plt.ylim(ymin,ymax)
    plt.xlabel('xi')
    plt.title('Mínimos cuadrados')
    
    plt.show()
    
  • Medir el tiempo de ejecución en Python

    Con varios métodos para resolver un problema, una media de eficiencia es el tiempo de ejecución de un algoritmo.

    dibujo cronómetro y banderaPara determinar el tiempo de ejecución de un algoritmo, se toman las lecturas de tiempo antes y después del bloque, calculando el tiempo de ocupación de las diferencias entre tiempos.

    La librería "time" permite obtener lectura del reloj del computador.

    Como un algoritmo de prueba se usa la suma de los m primeros números enteros.

    # tiempos de ejecuci n de algoritmos
    import time
    
    # INGRESO
    m = 100
    
    # PROCEDIMIENTO
    t1 = time.clock()
    
    # Ejecuta Operaciones del algoritmo
    suma = 0
    for i in range(0,m,1):
        suma = suma + i
    
    t2 = time.clock()
    # Tiempo para ejecutar operaciones
    ocupado = t2-t1
    
    # SALIDA
    print('tiempos:')
    print(t1, t2)
    print('tiempo de ejecuci n:')
    print(ocupado)
    

    con lo que se obtienen los siguientes resultados:

    tiempos:
    8.210955878428587e-07 5.46028565915501e-05
    tiempo de ejecución:
    5.378176100370724e-05
    >>> 
    
  • Sympy - Fórmulas y funciones simbólicas

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]

    SymPy es una librería usada para manejar de forma algebraica las expresiones matemáticas. Se requiere definir el símbolo que representa la variable en la expresión, por ejemplo la letra 'x'. Si la librería Sympy no está disponible o muestra un error la puede revisar las instrucciones del enlace instalar con pip.

    Una formula o función matemática f(x) descrita como fx se puede derivar, integrar, simplificar sus términos. También se puede construir una expresión matemática, por ejemplo desarrollar los términos para un polinomio como Taylor.

    Referencia: https://www.sympy.org/es/

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    1. Expresiones con Sympy

    Una función f(x) como un polinomio

    f(x) = (1-x)^5 + 5 x ((1-x)^4) - 0.4

    Se pueden manejar de forma algebraica con Sympy al declarar la variable independiente 'x' como un símbolo. La expresión se asigna a fx

    import sympy as sym
    x = sym.Symbol('x')
    fx = (1-x)**5 + 5*x*((1-x)**4) - 0.4
    

    Se simplifican o se expande los términos del polinomio con solo una instrucción,

    fx = fx.expand()
    print(fx)
    

    mostrando el siguiente resultado:

    4*x**5 - 15*x**4 + 20*x**3 - 10*x**2 + 0.6
    

    o una presentación diferente con sym.pprint():

    >>> sym.pprint(fx)
               4          5      
    5*x*(1 - x)  + (1 - x)  - 0.4
    >>> 
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/printing.html#setting-up-pretty-printing

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    2. Operaciones, combinar otras expresiones

    A una función simbólica fx se pueden añadir más términos con el mismo símbolo

    >>> import sympy as sym
    >>> sym.Symbol('x')
    >>> fx = sym.cos(x)
    >>> gx = fx + x
    >>> gx
    x + cos(x)
    >>> gx = gx - 2
    >>> gx
    x + cos(x) - 2
    

    Por lo que las funciones simbólicas son útiles cuando se construyen expresiones, como en el caso de series de Taylor descritas en la Unidad01

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    3. Evaluar una expresión f(x) con Sympy

    Las expresiones simbólicas se pueden evaluar en un punto, ejemplo x0=0.1 usando la instrucción fx.subs(x,x0), que sustituye el símbolo de la variable x con el valor definido para x0.

    >>> import sympy as sym
    >>> sym.Symbol('x')
    >>> fx = sym.cos(x)
    >>> fx.subs(x,0.1)
    0.995004165278026
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/basic_operations.html#substitution

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    4. Conversión de forma simbólica a forma numérica Lambda

    Para evaluar varios puntos en la expresión  en forma numérica para usar 'x' con valores en un vector o matriz como un arreglo, se convierte a la forma numérica Lambda con la instrucción sym.lambdify(x,fx)

    >>> >>> import sympy as sym
    >>> sym.Symbol('x')
    >>> fx = sym.cos(x)
    &>>> f_x = sym.lambdify(x,fx)
    >>> f_x(0.1)
    0.9950041652780258
    >>> f_x([0.1,0.3,0.5])
    array([0.99500417, 0.95533649, 0.87758256])
    >>> 
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/basic_operations.html#lambdify

    Ejemplo: Polinomio de Taylor – Ejemplos con Sympy-Python

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    5. expresiones de Funciones matemáticas

    Las Funciones matemáticas usadas en otras librerías como Numpy tienen también su representación simbólica en Sympy. Ejemplo cos(x):

    f(x) = \cos (x)
    >>> import sympy as sym
    >>> x = sym.Symbol('x')
    >>> fx = sym.cos(x)
    

    Realice pruebas con otras funciones: sin(x), exp(x), log(x), log10(x), Heaviside(x)

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    6. Derivadas con Sympy

    Las expresiones de la derivada se obtienen con la expresión fx.diff(x,k), indicando la k-ésima  derivada de la expresión.

    >>> import sympy as sym
    >>> x = sym.Symbol('x')
    >>> fx = sym.cos(x)
    >>> x
    x
    >>> fx
    cos(x)
    >>> fx.diff(x,1)
    -sin(x)
    >>> fx.diff(x,2)
    -cos(x)
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/calculus.html#derivatives

    Ejemplos: Sistemas LTI CT – Respuesta a entrada cero con Sympy-Python

    Derivadas como expresión de Sympy sin evaluar

    Cuando se requiere expresar tan solo la operación de derivadas, que será luego usada o reemplazada con otra expresión, se usa la derivada sin evaluar. Ejemplo:

    y = \frac{d}{dx}f(x)

    Para más adelante definir f(x).

    En Sympy, la expresión de y se realiza indicando que f será una variable

    x = sym.Symbol('x', real=True)
    f = sym.Symbol('f', real=True)
    y = sym.diff(f,x, evaluate=False) # derivada sin evaluar
    g = sym.cos(x) + x**2
    yg = y.subs(f,g).doit() # sustituye f con g y evalua .doit()
    
    >>> y
    Derivative(f, x)
    >>> g
    x**2 + cos(x)
    >>> yg
    2*x - sin(x)

    Ejemplos:

    Polinomio de Taylor – Ejemplos con Sympy-Python

    EDO lineal – solución complementaria y particular con Sympy-Python

    EDO lineal – ecuaciones auxiliar, general y complementaria con Sympy-Python

    EDP Parabólica – analítico con Sympy-Python

    EDP Elípticas – analítico con Sympy-Python

    Evaluación de las expresiones Sympy

    Se define la expresión 'derivada' y se la usa con la instrucción fx.subs(x,valor)

    >>> fx.subs(x,0)
    1
    >>> fx.subs(x,1/3)
    0.944956946314738
    >>> derivada = fx.diff(x,1)
    >>> derivada
    -sin(x)
    >>> derivada.subs(x,1/3)
    -0.327194696796152
    >>> 
    

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    7. Integrales con Sympy

    Sympy incorpora la operación del integral en sus expresiones, que pueden ser integrales definidas en un intervalo o expresiones sin evaluar.

    >>> import sympy as sym
    >>> t = sym.Symbol('t',real=True)
    >>> x = 10*sym.exp(-3*t)
    >>> x
    10*exp(-3*t)
    >>> y = sym.integrate(x,(t,0,10))
    >>> y
    10/3 - 10*exp(-30)/3
    >>> y = sym.integrate(x,(t,0,sym.oo))
    >>> y
    10/3
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/calculus.html#integrals

    Ejemplos:
    Transformada de Laplace – Integral con Sympy-Python

    Series de Fourier de señales periódicas con Sympy-Python

    Integral de Convolución entre x(t) y h(t) causal con Sympy-Python

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]

  • Numpy - Vectores y matrices

    Resumen de algunas funciones para vectores y matrices como arreglos en la librería Numpy de Python, usadas en el curso para facilitan el cálculo numérico. El orden de las instrucciones es el que aparece en las entradas del blog.

    Lo primero es hacer el llamado a las librerías con el alias 'np'

    import numpy as np
    

    Vectores

    puntos espaciados entre [a,b] - np.linspace()

    para obtener puntos en el rango [a,b] para una cantidad de tramos. Se obtienen tramos+1 puntos .

    a = 1
    b = 2
    tramos = 4
    x = np.linspace(a,b,tramos+1)
    
    >>> x
    array([ 1.  ,  1.25,  1.5 ,  1.75,  2.  ])
    

    rango de valores en vector - np.arange(a,b,dt)

    crea un vector con valores en el rango [a,b) y espaciados dt.

    t=np.arange(0,10,2)
    >>> t
    array([0, 2, 4, 6, 8])

    transponer vector - np.transpose()

    >>> fila = np.array([1,2,3])
    >>> columna = np.transpose([fila])
    >>> columna
    array([[1],
           [2],
           [3]])
    

    indices para ordenar - np.argsort()

    para ordenar por una columna específica:
    - se obtiene un vector con la columna que se requiere ordenar: tabla[:,0]
    - con el vector se determinan los índices para ordenar la tabla por filas: np.argsort()
    - se aplica los índices a una copia de la matriz: tabla[orden]

    tabla = np.array([[5,3],
                      [4,2],
                      [3,1],
                      [2,4],
                      [1,5]])
    
    # ordenar por primera columna
    referencia = tabla[:,0]
    orden = np.argsort(referencia)
    ordenada = tabla[orden]
    
    print(orden)
    print(ordenada)
    

    resultado:

    [4 3 2 1 0]
    [[1 5]
     [2 4]
     [3 1]
     [4 2]
     [5 3]]
    
    # ordenar por segunda columna
    referencia = tabla[:,1]
    orden = np.argsort(referencia)
    ordenada = tabla[orden]
    print(orden)
    print(ordenada)
    

    resultado:

    [2 1 0 3 4]
    [[3 1]
     [4 2]
     [5 3]
     [2 4]
     [1 5]]
    

    Matrices en Numpy

    concatenar- np.concatenate((A,b), axis=1)

    concatenar usa el parámetro axis: 0 para filas, 1 para columnas.

    import numpy as np
    cantidad = np.array([[4,2,5],
                         [2,5,8],
                         [5,4,3]])
    
    pagado = np.array([[60.70],
                       [92.90],
                       [56.30]])
    
    matriz = np.concatenate((cantidad,pagado),axis=1)
    
    >>> matriz
    array([[  4. ,   2. ,   5. ,  60.7],
           [  2. ,   5. ,   8. ,  92.9],
           [  5. ,   4. ,   3. ,  56.3]])
    

    resolver sistema de ecuaciones - np.linalg.solve(A,B)

    Resuelve el sistema de ecuaciones dado por una matriz A y un vector B. siendo, por ejemplo:

     0 =  c1 +  c2
    -5 = -c1 - 2c2
    
    c1 = -5 
    c2 = 5
    >>> A = [[ 1, 1],
    	 [-1,-2]]
    >>> B = [0,-5]
    >>> np.linalg.solve(A,B)
    array([-5.,  5.])
    >>>

    Otras instrucciones

    np.pi constante con valor π

    >>> np.pi
    3.141592653589793
    np.sin(t)
    np.cos(t)
    función trigonométrica en radianes. La variable t puede ser un escalar o un arreglo.

    >>> t=0.65
    >>> np.sin(0.65)
    0.60518640573603955
    >>> t=[0, 0.3, 0.6]
    >>> np.sin(t)
    array([ 0. , 0.29552021, 0.56464247])
    np.abs() obtiene el valor absoluto de un número. En el caso de un número complejo obtiene la parte real.
    np.real(complejo)
    np.imag(complejo)
    obtiene la parte real de los números complejos en un vector. Se aplica lo mismo para la parte imaginaria del número complejo.
    complex(a,b) crea el número complejo a partir de los valores de a y b.
    a=2
    b=3
    el resultado es: 2+3j
    np.piecewise(t, t>=donde, [1,0]) función que crea a partir de t, los valores de la condición t>=donde, ubicando los valores de 1, para otro caso es 0. Usada en la función escalón.
    np.roots([a,b,c]) obtiene las raíces del polinomio:
    ax2+bx+c
    x2 + 3 x + 2 = (x+1)(x+2)

    >>> np.roots([1,3,2])
    array([-2., -1.])
  • Funciones lambda vs def-return

    Las dos formas de escritura de funciones matemáticas básicamente hacen lo mismo. Sin embargo, cuando la función se define por tramos, la forma def-return se convierte en la más usada.

    Se adjunta un ejemplo:

    Funciones Lambda

    Permite describir funciones sencillas de una sola línea, que no está conformada por intervalos.

    f(x) = x \cos(x) - 2 x^2 + 3 x -1
    import numpy as np
    fx = lambda x: x*np.cos(x) - 2*(x**2) + 3*x -1
    
    >>> fx(0.2)
    -0.28398668443175157
    

    Funciones def-return

    Permiten describir en detalle lo que sucede con una función matemática por intervalos. Tiene la ventaja de permitir definir valores por intervalos, puntos discontínuos, etc.

    f(x)= \begin{cases} x \cos(x) - 2 x^2 + 3 x -1, & x>0 \\1, & x \leq 0 \end{cases}
    import numpy as np
    def fx(x):
        if x>0:
            fi = x*np.cos(x) - 2*(x**2) + 3*x -1
        if x<=0:
            fi = 1 
        return(fi)
    
    >>> fx(0.2)
     -0.28398668443175157
    
  • 7.3 EDP hiperbólicas con Python

    [ concepto ] [ analítico ] [ algoritmo ]
    ..


    1. Concepto y ejercicio

    Referencia:  Chapra PT8.1 p860,  Rodríguez 10.4 p435

    Las Ecuaciones Diferenciales Parciales tipo hiperbólicas semejantes a la mostrada, para un ejemplo en particular, representa la ecuación de onda de una dimensión u[x,t], que representa el desplazamiento vertical de una cuerda.

    \frac{\partial ^2 u}{\partial t^2}=c^2\frac{\partial ^2 u}{\partial x^2}

    EDP Cuerda 01

    Los extremos de la cuerda de longitud unitaria están sujetos y referenciados a una posición 0 a la izquierda y 1 a la derecha.

    u[x,t] , 0<x<1, t≥0
    u(0,t) = 0 , t≥0
    u(1,t) = 0 , t≥0

    Al inicio, la cuerda está estirada por su punto central:

    u(x,0) = \begin{cases} -0.5x &, 0\lt x\leq 0.5 \\ 0.5(x-1) &, 0.5\lt x \lt 1 \end{cases}

    Se suelta la cuerda, con velocidad cero desde la posición inicial:

    \frac{\partial u(x,0)}{\partial t} = 0

    [ concepto ] [ analítico ] [ algoritmo ]
    ..


    2. Desarrollo analítico

    La solución se realiza de forma semejante al procedimiento para EDP parabólicas y elípticas. Se cambia a la forma discreta  usando diferencias finitas divididas:

    \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta t)^2} =c^2 \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}

    El error es del orden O(\Delta x)^2 + O(\Delta t)^2.
    se reagrupa de la forma:

    u_{i,j+1}-2u_{i,j}+u_{i,j-1} = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \big( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \big)

    El término constante, por facilidad se simplifica con el valor de 1

    \lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} =1

    si c = 2 y Δx = 0.2, se deduce que Δt = 0.1

    que al sustituir en la ecuación, se simplifica anulando el término u[i,j]:

    u_{i,j+1}+u_{i,j-1} = u_{i+1,j}+u_{i-1,j}

    EDP Cuerda 02

    en los que intervienen solo los puntos extremos. Despejando el punto superior del rombo:

    u_{i,j+1} = u_{i+1,j}-u_{i,j-1}+u_{i-1,j}

    Convergencia:

    \lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \leq 1

    para simplificar aún más el problema, se usa la condición velocidad inicial de la cuerda igual a cero

    \frac{\delta u_{i,0}}{\delta t}=\frac{u_{i,1}-u_{i,-1}}{2\Delta t} = 0 u_{i,-1}=u_{i,1}

    que se usa para cuando el tiempo es cero, permite calcular los puntos para t[1]:

    j=0

    u_{i,1} = u_{i+1,0}-u_{i,-1}+u_{i-1,0} u_{i,1} = u_{i+1,0}-u_{i,1}+u_{i-1,0} 2 u_{i,1} = u_{i+1,0}+u_{i-1,0} u_{i,1} = \frac{u_{i+1,0}+u_{i-1,0}}{2}

    Aplicando solo cuando j = 0

    EDP Cuerda 03

    que al ponerlos en la malla se logra un sistema explícito para cada u[i,j], con lo que se puede resolver con un algoritmo:

    EDP hiperbolica 01

    [ concepto ] [ analítico ] [ algoritmo ]
    ..


    3. Algoritmo con Python

    # Ecuaciones Diferenciales Parciales
    # Hiperbólica. Método explicito
    import numpy as np
    
    def cuerdainicio(xi):
        n = len(xi)
        y = np.zeros(n,dtype=float)
        for i in range(0,n,1):
            if (xi[i]<=0.5):
                y[i]=-0.5*xi[i]
            else:
                y[i]=0.5*(xi[i]-1)
        return(y)
    
    # INGRESO
    x0 = 0 # Longitud de cuerda
    xn = 1
    y0 = 0 # Puntos de amarre
    yn = 0
    t0 = 0 # tiempo inicial
    c = 2  # constante EDP
    # discretiza
    tramosx = 16
    tramost = 32
    dx = (xn-x0)/tramosx 
    dt = dx/c
    
    # PROCEDIMIENTO
    xi = np.arange(x0,xn+dx/2,dx)
    tj = np.arange(0,tramost*dt+dt/2,dt)
    n = len(xi)
    m = len(tj)
    
    u = np.zeros(shape=(n,m),dtype=float)
    u[:,0] = cuerdainicio(xi)
    u[0,:] = y0
    u[n-1,:] = yn
    # Aplicando condición inicial
    j = 0
    for i in range(1,n-1,1):
        u[i,j+1] = (u[i+1,j]+u[i-1,j])/2
    # Para los otros puntos
    for j in range(1,m-1,1):
        for i in range(1,n-1,1):
            u[i,j+1] = u[i+1,j]-u[i,j-1]+u[i-1,j]
    
    # SALIDA
    np.set_printoptions(precision=2)
    print('xi =')
    print(xi)
    print('tj =')
    print(tj)
    print('matriz u =')
    print(u)
    

    con lo que se obtienen los resultados numéricos, que para mejor interpretación se presentan en una gráfica estática y otra animada.

    # GRAFICA
    import matplotlib.pyplot as plt
    
    for j in range(0,m,1):
        y = u[:,j]
        plt.plot(xi,y)
    
    plt.title('EDP hiperbólica')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.show()
    
    # **** GRÁFICO CON ANIMACION ***********
    import matplotlib.animation as animation
    
    # Inicializa parametros de trama/foto
    retardo = 70   # milisegundos entre tramas
    tramas  = m
    maximoy = np.max(np.abs(u))
    figura, ejes = plt.subplots()
    plt.xlim([x0,xn])
    plt.ylim([-maximoy,maximoy])
    
    # lineas de función y polinomio en gráfico
    linea_poli, = ejes.plot(xi,u[:,0], '-')
    plt.axhline(0, color='k')  # Eje en x=0
    
    plt.title('EDP hiperbólica')
    # plt.legend()
    # txt_x = (x0+xn)/2
    txt_y = maximoy*(1-0.1)
    texto = ejes.text(x0,txt_y,'tiempo:',
                      horizontalalignment='left')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    
    # Nueva Trama
    def unatrama(i,xi,u):
        # actualiza cada linea
        linea_poli.set_ydata(u[:,i])
        linea_poli.set_xdata(xi)
        linea_poli.set_label('tiempo linea: '+str(i))
        texto.set_text('tiempo['+str(i)+']')
        # color de la línea
        if (i<=9):
            lineacolor = 'C'+str(i)
        else:
            numcolor = i%10
            lineacolor = 'C'+str(numcolor)
        linea_poli.set_color(lineacolor)
        return linea_poli, texto
    
    # Limpia Trama anterior
    def limpiatrama():
        linea_poli.set_ydata(np.ma.array(xi, mask=True))
        linea_poli.set_label('')
        texto.set_text('')
        return linea_poli, texto
    
    # Trama contador
    i = np.arange(0,tramas,1)
    ani = animation.FuncAnimation(figura,
                                  unatrama,
                                  i ,
                                  fargs=(xi,u),
                                  init_func=limpiatrama,
                                  interval=retardo,
                                  blit=True)
    # Graba Archivo video y GIFAnimado :
    # ani.save('EDP_hiperbólica.mp4')
    ani.save('EDP_hiperbolica.gif', writer='imagemagick')
    plt.draw()
    plt.show()
    

    Una vez realizado el algoritmo, se pueden cambiar las condiciones iniciales de la cuerda y se observan los resultados.

    Se recomienda realizar otros ejercicios de la sección de evaluaciones para EDP Hiperbólicas y observar los resultados con el algoritmo modificado.

    [ concepto ] [ analítico ] [ algoritmo ]

  • 7.2.2 EDP Elípticas método implícito con Python

    EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]
    ..


    1. EDP Elípticas: Método Implícito – Desarrollo Analítico

    con el resultado desarrollado en EDP elípticas para:

    \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

    y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1

    se puede plantear que:

    u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

    EDP Elipticas Iterativo

    con lo que para el método implícito, se plantea un sistema de ecuaciones para determinar los valores en cada punto desconocido.

    j=1, i =1

    u_{2,1}-4u_{1,1}+u_{0,1} + u_{1,2} +u_{1,0} = 0 u_{2,1}-4u_{1,1}+Ta + u_{1,2} +Tc= 0 -4u_{1,1}+u_{2,1}+u_{1,2} = -(Tc+Ta)

    j=1, i =2

    u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +u_{2,0} = 0 u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +Tc = 0 u_{1,1}-4u_{2,1}+u_{3,1}+ u_{2,2}= -Tc

    j=1, i=3

    u_{4,1}-4u_{3,1}+u_{2,1} + u_{3,2} +u_{3,0} = 0 Tb-4u_{3,1}+u_{2,1} + u_{3,2} +Tc = 0 u_{2,1} -4u_{3,1} + u_{3,2} = -(Tc+Tb)

    j=2, i=1

    u_{2,2}-4u_{1,2}+u_{0,2} + u_{1,3} +u_{1,1} = 0 u_{2,2}-4u_{1,2}+Ta + u_{1,3} +u_{1,1} = 0 -4u_{1,2}+u_{2,2}+u_{1,1}+u_{1,3} = -Ta

    j = 2, i = 2

    u_{1,2}-4u_{2,2}+u_{3,2} + u_{2,3} +u_{2,1} = 0

    j = 2, i = 3

    u_{4,2}-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 Tb-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 u_{2,2} -4u_{3,2}+ u_{3,3} +u_{3,1} = -Tb

    j=3, i = 1

    u_{2,3}-4u_{1,3}+u_{0,3} + u_{1,4} +u_{1,2} = 0 u_{2,3}-4u_{1,3}+Ta + Td +u_{1,2} = 0 -4u_{1,3}+u_{2,3}+u_{1,2} = -(Td+Ta)

    j=3, i = 2

    u_{3,3}-4u_{2,3}+u_{1,3} + u_{2,4} +u_{2,2} = 0 u_{3,3}-4u_{2,3}+u_{1,3} + Td +u_{2,2} = 0 +u_{1,3} -4u_{2,3}+u_{3,3} +u_{2,2} = -Td

    j=3, i=3

    u_{4,3}-4u_{3,3}+u_{2,3} + u_{3,4} +u_{3,2} = 0 Tb-4u_{3,3}+u_{2,3} + Td +u_{3,2} = 0 u_{2,3}-4u_{3,3}+u_{3,2} = -(Td+Tb)

    con las ecuaciones se arma una matriz:

    A = np.array([
        [-4, 1, 0, 1, 0, 0, 0, 0, 0],
        [ 1,-4, 1, 0, 1, 0, 0, 0, 0],
        [ 0, 1,-4, 0, 0, 1, 0, 0, 0],
        [ 1, 0, 0,-4, 1, 0, 1, 0, 0],
        [ 0, 1, 0, 1,-4, 1, 0, 1, 0],
        [ 0, 0, 1, 0, 1,-4, 0, 0, 1],
        [ 0, 0, 0, 1, 0, 0,-4, 1, 0],
        [ 0, 0, 0, 0, 1, 0, 1,-4, 1],
        [ 0, 0, 0, 0, 0, 1, 0, 1,-4],
        ])
    B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb),
                      -Ta,   0,    -Tb,
                  -(Td+Ta),-Td,-(Td+Tb)])
    

    que al resolver el sistema de ecuaciones se obtiene:

    >>> Xu
    array([ 56.43,  55.71,  56.43,  60.  ,  60.  ,  60.  ,  63.57,  64.29,
            63.57])
    

    ingresando los resultados a la matriz u:

    xi=
    [ 0.   0.5  1.   1.5  2. ]
    yj=
    [ 0.    0.38  0.75  1.12  1.5 ]
    matriz u
    [[ 60.    60.    60.    60.    60.  ]
     [ 50.    56.43  60.    63.57  70.  ]
     [ 50.    55.71  60.    64.29  70.  ]
     [ 50.    56.43  60.    63.57  70.  ]
     [ 60.    60.    60.    60.    60.  ]]
    >>>
    

    EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]

    ..


    2. EDP Elípticas: Método Implícito – Desarrollo con algoritmo

    Instrucciones en Python

    # Ecuaciones Diferenciales Parciales
    # Elipticas. Método implícito
    import numpy as np
    
    # INGRESO
    # Condiciones iniciales en los bordes
    Ta = 60
    Tb = 60
    Tc = 50
    Td = 70
    # dimensiones de la placa
    x0 = 0
    xn = 2
    y0 = 0
    yn = 1.5
    # discretiza, supone dx=dy
    tramosx = 4
    tramosy = 4
    dx = (xn-x0)/tramosx 
    dy = (yn-y0)/tramosy 
    maxitera = 100
    tolera = 0.0001
    
    A = np.array([
        [-4, 1, 0, 1, 0, 0, 0, 0, 0],
        [ 1,-4, 1, 0, 1, 0, 0, 0, 0],
        [ 0, 1,-4, 0, 0, 1, 0, 0, 0],
        [ 1, 0, 0,-4, 1, 0, 1, 0, 0],
        [ 0, 1, 0, 1,-4, 1, 0, 1, 0],
        [ 0, 0, 1, 0, 1,-4, 0, 0, 1],
        [ 0, 0, 0, 1, 0, 0,-4, 1, 0],
        [ 0, 0, 0, 0, 1, 0, 1,-4, 1],
        [ 0, 0, 0, 0, 0, 1, 0, 1,-4],
        ])
    B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb),
                  -Ta,0,-Tb,
                  -(Td+Ta),-Td,-(Td+Tb)])
    
    
    # PROCEDIMIENTO
    # Resuelve sistema ecuaciones
    Xu = np.linalg.solve(A,B)
    [nx,mx] = np.shape(A)
    
    xi = np.arange(x0,xn+dx/2,dx)
    yj = np.arange(y0,yn+dy/2,dy)
    n = len(xi)
    m = len(yj)
    
    u = np.zeros(shape=(n,m),dtype=float)
    u[:,0]   = Tc
    u[:,m-1] = Td
    u[0,:]   = Ta
    u[n-1,:] = Tb
    u[1:1+3,1] = Xu[0:0+3]
    u[1:1+3,2] = Xu[3:3+3]
    u[1:1+3,3] = Xu[6:6+3]
    
    # SALIDA
    np.set_printoptions(precision=2)
    print('xi=')
    print(xi)
    print('yj=')
    print(yj)
    print('matriz u')
    print(u)
    

    La gráfica de resultados se obtiene de forma semejante al ejercicio con método iterativo.

    Se podría estandarizar un poco más el proceso para que sea realizado por el algoritmo y sea más sencillo generar la matriz con más puntos. Tarea.

  • 7.2.1 EDP Elípticas método iterativo con Python

    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]

    ..


    1. Ejercicio

    Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

    Siguiendo el tema anterior en EDP Elípticas, para resolver la parte numérica asuma como:

    Valores de frontera: Ta = 60 , Tb = 25 , Tc = 50, Td = 70
    Longitud en:  x0 = 0, xn = 2 , y0 = 0, yn = 1.5
    Tamaño de paso dx = 0.25 dy = dx
    iteramax = 100 tolera = 0.0001

    Para la ecuación diferencial parcial Elíptica

    \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0


    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]
    ..


    2. EDP Elípticas: Método iterativo - Desarrollo Analítico

    A partir del resultado desarrollado en EDP elípticas:

    \frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda u_{i+1,j} +(-2\lambda-2) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1} +u_{i,j-1}=0

    y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1

    se puede plantear que:

    u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

    que reordenando para un punto central desconocido se convierte a:

    u_{i,j} = \frac{1}{4} \big[ u_{i+1,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} \big]

    con lo que se interpreta que cada punto central es el resultado del promedio de los puntos alrededor del rombo formado en la malla.

    El cálculo numérico se puede realizar de forma iterativa haciendo varias pasadas en la matriz, promediando cada punto. Para revisar las iteraciones se controla la convergencia junto con un máximo de iteraciones.

    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]

    ..


    3. Algoritmo en Python

    Para el ejercicio dado, se presenta el resultado para tres iteraciones y el resultado final con un gráfico en 3D:

    u inicial:
    [[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
     [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
     [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
     [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
     [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
     [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
     [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
    itera: 0   ;  u:
    [[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
     [60.   53.12 51.41 50.98 50.87 50.84 50.84 44.27 25.  ]
     [60.   53.91 51.95 51.36 51.18 51.13 51.12 42.91 25.  ]
     [60.   54.1  52.14 51.5  51.3  51.23 51.21 42.59 25.  ]
     [60.   54.15 52.2  51.55 51.34 51.27 51.24 42.52 25.  ]
     [60.   58.85 58.07 57.72 57.58 57.52 57.5  48.76 25.  ]
     [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
    errado_u:  8.728094100952148
    itera: 1   ;  u:
    [[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
     [60.   53.83 51.69 50.98 50.75 50.68 49.02 41.73 25.  ]
     [60.   54.97 52.54 51.55 51.18 51.05 48.55 39.47 25.  ]
     [60.   55.31 52.89 51.82 51.39 51.23 48.4  38.85 25.  ]
     [60.   56.59 54.78 53.91 53.54 53.38 50.45 40.76 25.  ]
     [60.   61.17 60.91 60.6  60.42 60.33 57.38 48.29 25.  ]
     [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
    errado_u:  3.7443900108337402
    itera: 2   ;  u:
    [[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
     [60.   54.17 51.92 51.06 50.73 50.2  47.62 40.52 25.  ]
     [60.   55.5  52.97 51.76 51.23 50.3  46.45 37.7  25.  ]
     [60.   56.25 53.95 52.75 52.19 51.07 46.71 37.54 25.  ]
     [60.   58.05 56.71 55.9  55.47 54.33 49.8  40.16 25.  ]
     [60.   62.24 62.39 62.18 61.99 60.93 57.25 48.1  25.  ]
     [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
    errado_u:  2.0990657806396484
    ... continua
    Método Iterativo EDP Elíptica
    iteraciones: 41  converge =  True
    xi: [0.   0.25 0.5  0.75 1.   1.25 1.5  1.75 2.  ]
    yj: [0.   0.25 0.5  0.75 1.   1.25 1.5 ]
    Tabla de resultados en malla EDP Elíptica
    j, U[i,j]
    6 [70. 70. 70. 70. 70. 70. 70. 70. 70.]
    5 [60.   64.02 64.97 64.71 63.62 61.44 57.16 47.96 25.  ]
    4 [60.   61.1  61.14 60.25 58.35 54.98 49.23 39.67 25.  ]
    3 [60.   59.23 58.25 56.8  54.53 50.89 45.13 36.48 25.  ]
    2 [60.   57.56 55.82 54.19 52.09 48.92 43.91 36.14 25.  ]
    1 [60.   55.21 53.27 52.05 50.73 48.78 45.46 39.15 25.  ]
    0 [50. 50. 50. 50. 50. 50. 50. 50. 50.]
    >>>
    

    cuyos valores se interpretan mejor en una gráfica, en este caso 3D:

    EDP Elipticas Iterativo

    Instrucciones en Python

    # Ecuaciones Diferenciales Parciales
    # Elipticas. Método iterativo
    # d2u/dx2  + du/dt = f(x,y)
    import numpy as np
    
    # INGRESO
    fxy = lambda x,y: 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson
    # Condiciones iniciales en los bordes, valores de frontera
    Ta = 60     # izquierda de la placa
    Tb = 25     # derecha de la placa
    #Tc = 50    # inferior 
    fxc = lambda x: 50+0*x #  f(x) inicial inferior
    # Td = 70   # superior
    fxd = lambda x: 70+0*x #  f(x) inicial superior
    
    # dimensiones de la placa
    x0 = 0    # longitud en x
    xn = 2
    y0 = 0    # longitud en y
    yn = 1.5
    # discretiza, supone dx=dy
    dx = 0.25  # Tamaño de paso
    dy = dx
    
    iteramax = 100  # maximo de iteraciones
    tolera = 0.0001
    verdigitos = 2   # decimales a mostrar en tabla de resultados
    vertabla = True  # ver iteraciones
    
    # coeficientes de U[x,t]. factores P,Q,R
    lamb = (dy**2)/(dx**2)
    P = lamb       # izquierda P*U[i-1,j]
    Q = -2*lamb-2  # centro Q*U[i,j]
    R = lamb       # derecha R*U[i+1,j]
    
    # PROCEDIMIENTO
    # iteraciones en xi,yj ancho,profundidad
    xi = np.arange(x0,xn+dx/2,dx)
    yj = np.arange(y0,yn+dy/2,dy)
    n = len(xi)
    m = len(yj)
    
    # Matriz u[xi,yj], tabla de resultados
    u = np.zeros(shape=(n,m),dtype=float)
    
    # llena u con valores en fronteras
    u[0,:]   = Ta   # izquierda
    u[n-1,:] = Tb   # derecha
    fic = fxc(xi)   # Tc inferior
    u[:,0]   = fic
    fid = fxd(xi)   # Td superior
    u[:,m-1] = fid  
    # estado inicial dentro de la placa
    # promedio = (Ta+Tb+Tc+Td)/4
    promedio = (Ta+Tb+np.max(fic)+np.max(fid))*(1/-Q)
    u[1:n-1,1:m-1] = promedio
    
    np.set_printoptions(precision=verdigitos)
    if vertabla == True:
        print('u inicial:')
        print(np.transpose(u))
    
    # iterar puntos interiores
    U_3D = [np.copy(u)]
    U_error = ['--']
    itera = 0
    converge = False
    while (itera<=iteramax and converge==False):
        itera = itera +1
        Uantes = np.copy(u)
        for i in range(1,n-1):
            for j in range(1,m-1):
                u[i,j] = P*u[i-1,j]+R*u[i+1,j]+u[i,j-1]+u[i,j+1]
                u[i,j] = u[i,j] - fxy(xi[i],yj[j])*(dy**2)
                u[i,j] = (1/-Q)*u[i,j]
        diferencia = u - Uantes
        errado_u = np.max(np.abs(diferencia))
        if (errado_u<tolera):
            converge=True
        U_3D.append(np.copy(u))
        U_error.append(np.around(errado_u,verdigitos))
        
        if (vertabla==True) and (itera<=3):
            np.set_printoptions(precision=verdigitos)
            print('itera:',itera-1,'  ; ','u:')
            print(np.transpose(u))
            print('errado_u: ', errado_u)
        if (vertabla==True) and (itera==4):
            print('... continua')
            
    # SALIDA
    print('Método Iterativo EDP Elíptica')
    print('iteraciones:',itera,' converge = ', converge)
    print('errado_u: ', errado_u)
    print('xi:', xi)
    print('yj:', yj)
    print('Tabla de resultados en malla EDP Elíptica')
    print('j, U[i,j]')
    for j in range(m-1,-1,-1):
        print(j,u[:,j])
    

    La gráfica de resultados requiere ajuste de ejes, pues el índice de filas es el eje x, y las columnas es el eje y. La matriz y sus datos en la gráfica se obtiene como la transpuesta de u

    # GRAFICA 3D ------
    import matplotlib.pyplot as plt
    # Malla para cada eje X,Y
    Xi, Yi = np.meshgrid(xi, yj)
    U = np.transpose(u) # ajuste de índices fila es x
    
    fig_3D = plt.figure()
    graf_3D = fig_3D.add_subplot(111, projection='3d')
    graf_3D.plot_wireframe(Xi,Yi,U,
                           color ='blue')
    graf_3D.plot(Xi[1,0],Yi[1,0],U[1,0],'o',color ='orange')
    graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='red',
                 label='U[i,j]')
    graf_3D.plot(Xi[1,2],Yi[1,2],U[1,2],'o',color ='salmon')
    graf_3D.plot(Xi[0,1],Yi[0,1],U[0,1],'o',color ='green')
    graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')
    
    graf_3D.set_title('EDP elíptica')
    graf_3D.set_xlabel('x')
    graf_3D.set_ylabel('y')
    graf_3D.set_zlabel('U')
    graf_3D.legend()
    graf_3D.view_init(35, -45)
    plt.show()
    

    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]

     

  • 7.2 EDP Elípticas

    EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]
    ..


    1. EDP Elípticas

    Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

    Las Ecuaciones Diferenciales Parciales tipo elípticas semejantes a la mostrada:

    \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

    (ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)

    se interpreta como una placa metálica de dimensiones Lx y Ly, delgada con aislante que recubren las caras de la placa, y sometidas a condiciones en las fronteras:

    Lx = dimensión x de placa metálica
    Ly = dimensión y de placa metálica
    u[0,y]  = Ta
    u[Lx,y] = Tb
    u[x,0]  = Tc
    u[x,Ly] = Td

    Placa Metalica 01

    Para el planteamiento se usa una malla en la que cada nodo corresponden a los valores u[xi,yj]. Para simplificar la nomenclatura se usan los subíndices i para el eje de las x,  j para el eje t, quedando u[i,j].

    EDP Elipticas Iterativo

    vista superior, plano xy:

    Placa Metalica 02

    La ecuación se cambia a la forma discreta, usando diferencias divididas que se sustituyen en la ecuación, por ejemplo:

    \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta y)^2}=0

    Se agrupan los términos de los diferenciales:

    \frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda u_{i+1,j} - 2\lambda u_{i,j} + \lambda u_{i-1,j} + u_{i,j+1} - 2 u_{i,j} + u_{i,j-1} = 0 \lambda u_{i+1,j} + (-2\lambda-2) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1} +u_{i,j-1}=0

    Obteniendo así la solución numérica conceptual. La forma de resolver el problema determina el nombre del método a seguir.

    Una forma de simplificar la expresión, es hacer λ=1, es decir  \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1, se determina que los tamaño de paso Δx, Δy son iguales.

    u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

    EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]