Categoría: Unidades

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}

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

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}

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

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:

Algoritmo en 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
xn = 1 # Longitud de cuerda
y0 = 0
yn = 0 # Puntos de amarre
t0 = 0
c = 2
# discretiza
tramosx = 16
tramost = 32
dx = (xn-x0)/tramosx 
dt = dx/c

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
tj = np.arange(0,tramost*dt+dt,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.

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

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 ]

..

---

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

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,dx)
yj = np.arange(y0,yn+dy,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.

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

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

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:

Instrucciones en Python

# Ecuaciones Diferenciales Parciales
# Elipticas. Método iterativo
import numpy as np

# INGRESO
# Condiciones iniciales en los bordes
Ta = 60     # izquierda
Tb = 25     # derecha
#Tc = 50    # inferior 
fxc = lambda x: 50+0*x # función inicial inferior
# Td = 70   # superior
fxd = lambda x: 70+0*x # función 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
tolera = 0.0001
verdigitos = 2   # decimales a mostrar en tabla de resultados

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)

# Matriz u
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)
u[:,0]   = fic  # Tc inferior
fid = fxd(xi)
u[:,m-1] = fid  # Td superior
# valor inicial de iteración dentro de u
# promedio = (Ta+Tb+Tc+Td)/4
promedio = (Ta+Tb+np.max(fic)+np.max(fid))/4
u[1:n-1,1:m-1] = promedio
np.set_printoptions(precision=verdigitos)
print('u inicial:')
print(np.transpose(u))

# iterar puntos interiores
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] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1])/4
    diferencia = u - Uantes
    errado_u = np.max(np.abs(diferencia))
    if (errado_u<tolera):
        converge=True
    if itera<=3:
        np.set_printoptions(precision=verdigitos)
        print('itera:',itera-1,'  ; ','u:')
        print(np.transpose(u))
        print('errado_u: ', errado_u)
    if itera==4:
        print('... continua')
        
# SALIDA
print('Método Iterativo EDP Elíptica')
print('iteraciones:',itera,' converge = ', converge)
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 en 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 ]

 

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

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

Se discretiza la ecuación 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

con lo que se simplifican los valores como uno solo \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1 . Por facilidad de lo que se realiza se supone que lambda tiene valor de 1 o los delta son iguales.

u_{i+1,j}-4u_{i,j}+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.

---

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

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]

..

---

1. Ejercicio

Referencia:  Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodríguez 10.2.4 p415

Siguiendo el tema anterior en EDP Parabólicas, se tiene que:

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

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]

..

---

2. Desarrollo Analítico

En éste caso se usan diferencias finitas centradas y hacia atrás; la línea de referencia es t₁:

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

La selección de las diferencias divididas corresponden a los puntos que se quieren usar para el cálculo, se observan mejor en la gráfica de la malla.

Luego se sustituyen en la ecuación del problema, obteniendo:

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

De la gráfica se destaca que en la fórmula, dentro del triángulo solo hay DOS valores desconocidos, destacados por los punto en amarillo.
En la ecuación se representa por U[i,j] y U[i+1,j]. Por lo que será necesario crear un sistema de ecuaciones sobre toda la línea de tiempo t₁ para resolver el problema.

Despejando la ecuación, se agrupan términos constantes: λ = \frac{\Delta t}{K (\Delta x)^2} .

\lambda u_{i-1,j} + (-1-2\lambda) u_{i,j} + \lambda u_{i+1,j} = -u_{i,j-1}

Los parámetro P, Q y R se determinan de forma semejante al método explícito:

P = \lambda Q = -1-2\lambda R = \lambda Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j} = -u_{i,j-1}

Los valores en los extremos son conocidos, para los puntos intermedios  se crea un sistema de ecuaciones para luego usar la forma Ax=B y resolver los valores para cada u(x_i,t_j).

Por ejemplo con cuatro tramos entre extremos se tiene que:
indice de tiempo es 1 e índice de x es 1.

i=1,j=1

Pu_{0,1} + Qu_{1,1} + Ru_{2,1} = -u_{1,0}

i=2,j=1

Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -u_{2,0}

i=3,j=1

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

agrupando ecuaciones y sustituyendo valores conocidos:
\begin{cases} Qu_{1,1} + Ru_{2,1} + 0 &= -T_{0}-PT_{A}\\Pu_{1,1} + Qu_{2,1} + Ru_{3,1} &= -T_{0}\\0+Pu_{2,1}+Qu_{3,1}&=-T_{0}-RT_{B}\end{cases}

que genera la matriz a resolver:

\begin{bmatrix} Q && R && 0 && -T_{0}-PT_{A}\\P && Q && R && -T_{0}\\0 && P && Q && -T_{0}-RT_{B}\end{bmatrix}

Use alguno de los métodos de la unidad 3 para resolver el sistema y obtener los valores correspondientes.

Por la extensión de la solución es conveniente usar un algoritmo y convertir los pasos o partes pertinentes a funciones.

Tarea: Revisar y comparar con el método explícito.

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]

..

---

Algoritmos en Python

Para la solución con el método implícito, se obtiene el mismo resultado en la gráfica y tabla. Aunque el algoritmo es diferente.

algunos valores:

iteración j: 1
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-40. -25. -25. -25. -25. -25. -25. -25. -35.]
resultado en j: 1 
 [31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57]
iteración j: 2
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-46.01 -26.03 -25.18 -25.03 -25.01 -25.01 -25.08 -25.44 -37.57]
resultado en j: 2 
 [35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39]
iteración j: 3
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-50.25 -27.49 -25.55 -25.12 -25.03 -25.05 -25.24 -26.07 -39.39]
resultado en j: 3 
 [38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72]
EDP Parabólica - Método implícito 
lambda:  0.25
x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ] ... 1.0
Tabla de resultados en malla EDP Parabólica
j, U[:,: 10 ], primeras iteraciones
10 [60.   47.43 37.58 31.3  27.98 26.64 26.63 27.83 30.43 34.63 40.  ]
9 [60.   46.75 36.68 30.56 27.48 26.3  26.33 27.47 30.04 34.33 40.  ]
8 [60.   45.96 35.69 29.8  27.01 26.   26.06 27.12 29.6  33.99 40.  ]
7 [60.   45.01 34.6  29.02 26.57 25.73 25.81 26.76 29.13 33.58 40.  ]
6 [60.   43.87 33.39 28.24 26.16 25.51 25.58 26.41 28.6  33.09 40.  ]
5 [60.   42.45 32.06 27.48 25.8  25.32 25.39 26.07 28.03 32.48 40.  ]
4 [60.   40.67 30.61 26.75 25.51 25.18 25.24 25.76 27.41 31.71 40.  ]
3 [60.   38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72 40.  ]
2 [60.   35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39 40.  ]
1 [60.   31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
>>> 

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método implícito
# Referencia: Chapra 30.3 p.895 pdf.917
#       Rodriguez 10.2.5 p.417
import numpy as np

# INGRESO
# Valores de frontera
Ta = 60
Tb = 40
#T0 = 25 # estado inicial de barra
fx = lambda x: 25 + 0*x # función inicial para T0
a = 0  # longitud en x
b = 1

K = 4     # Constante K
dx = 0.1  # Tamaño de paso
dt = dx/10

n = 100 # iteraciones en tiempo
verdigitos = 2   # decimales a mostrar en tabla de resultados

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
m  = len(xi)
ultimox = m-1
ultimot = n-1
# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j = 0
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = -1 -2*lamb
R = lamb

# Calcula U para cada tiempo + dt
tj = np.arange(0,(n+1)*dt,dt)
j = 1
while not(j>=n):
    # Matriz de ecuaciones
    k = m-2
    A = np.zeros(shape=(k,k), dtype = float)
    B = np.zeros(k, dtype = float)
    for f in range(0,k,1):
        if (f>0):
            A[f,f-1]=P
        A[f,f] = Q
        if (f<(k-1)):
            A[f,f+1]=R
        B[f] = -u[f+1,j-1]
    B[0] = B[0]-P*u[0,j]
    B[k-1] = B[k-1]-R*u[m-1,j]
    # Resuelve sistema de ecuaciones
    C = np.linalg.solve(A, B)

    # copia resultados a u[i,j]
    for f in range(0,k,1):
        u[f+1,j] = C[f]

    # siguiente iteración
    j = j + 1

    # muestra 3 primeras iteraciones
    np.set_printoptions(precision=verdigitos)
    if j<(3+2): 
        print('iteración j:',j-1)
        print('A:\n',A)
        print('B:\n',B)
        print('resultado en j:',j-1,'\n',C)
        
# SALIDA
print('EDP Parabólica - Método implícito ')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:len(xi)],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',ultimox,'], primeras iteraciones')
for j in range(ultimox,-1,-1):
    print(j,u[:,j])

Para realizar la gráfica se aplica lo mismo que en el método explícito

# GRAFICA ------------
import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneas
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
plt.show()

Sin embargo para la gráfica en 3D se ajustan los puntos de referencia agregados por las diferencias finitas.

# GRAFICA en 3D
# vector de tiempo, simplificando en tramos
tramos = 10
salto = int(n/tramos)
if (salto == 0):
    salto = 1
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)

# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(tramos,m),dtype=float)
for k in range(0,tramos,1):
    U[k,:] = u[:,k*salto]
    tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)

fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U,
                       color ='blue',label='EDP Parabólica')
graf_3D.plot(Xi[2,0],Yi[2,0],U[2,0],'o',color ='orange')
graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')
graf_3D.plot(Xi[2,2],Yi[2,2],U[2,2],'o',color ='salmon')
graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='green')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

Queda por revisar la convergencia y estabilidad de la solución a partir de los O(h) de cada aproximación usada. Revisar los criterios en Chapra 30.2.1 p891, Burden 9Ed 12.2 p727, Rodríguez 10.2.2 p409 .

Tarea o proyecto: Realizar la comparación de tiempos de ejecución entre los métodos explícitos e implícitos. La parte de animación funciona igual en ambos métodos.  Los tiempos de ejecución se determinan usando las instrucciones descritas en el enlace: Tiempos de Ejecución en Python

---

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]

EDP Parabólica [ concepto ] [ ejercicio ]
Método explícito: [ Analítico ] [ Algoritmo ]

..

---

1. Ejercicio

Referencia:  Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodríguez 10.2 p406

Siguiendo el tema anterior en EDP Parabólicas, para resolver la parte numérica asuma

Valores de frontera: Ta = 60, Tb = 40, T0 = 25
Longitud en x a = 0, b = 1,  Constante K= 4
Tamaño de paso dx = 0.1, dt = dx/10

Para la ecuación diferencial parcial parabólica:

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

..

2. Método Explícito

Se usan las diferencias divididas, donde se requieren dos valores para la derivada de orden uno y tres valores para la derivada de orden dos:

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

La selección de las diferencias divididas corresponden a los puntos que se quieren usar para el cálculo, se observan mejor en la gráfica de la malla. La línea de referencia es el tiempo t₀

Luego se sustituyen en la ecuación del problema, obteniendo:

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

De la gráfica se destaca que en la fórmula solo hay un valor desconocido, destacado por el punto en amarillo dentro del triángulo. En la ecuación se representa por U[i,j+1].

Despejando la ecuación, se agrupan términos constantes:

λ = \frac{\Delta t}{K (\Delta x)^2}

quedando la ecuación, con los términos ordenados de izquierda a derecha como en la gráfica:

u_{i,j+1} = \lambda u_{i-1,j} +(1-2\lambda)u_{i,j}+\lambda u_{i+1,j}

Al resolver se encuentra que que cada valor en un punto amarillo se calcula como una suma ponderada de valores conocidos, por lo que el desarrollo se conoce como el método explícito.

La ponderación está determinada por los términos P, Q, y R.

\lambda = \frac{\Delta t}{K(\Delta x)^2} P = \lambda Q = 1-2\lambda R = \lambda u_{i ,j+1} = Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j}

Fórmulas que se desarrollan usando un algoritmo y considerando que al disminuir los valores de Δx y Δt la cantidad de operaciones aumenta.

Queda por revisar la convergencia y estabilidad de la solución a partir de los O(h) de cada aproximación usada.

Revisar los criterios en:  Chapra 30.2.1 p891 pdf915, Burden 9Ed 12.2 p727, Rodriguez 10.2.2 pdf409 .

\lambda \leq \frac{1}{2}

Cuando λ ≤ 1/2 se tiene como resultado una solución donde los errores no crecen, sino que oscilan.
Haciendo λ ≤ 1/4 asegura que la solución no oscilará.
También se sabe que con λ= 1/6 se tiende a minimizar los errores por truncamiento (véase Carnahan y cols., 1969).

Para resolver la parte numérica asuma que:

Valores de frontera: Ta = 60, Tb = 40, T0 = 25 
Longitud en x a = 0, b = 1 
Constante K= 4 
Tamaño de paso dx = 0.1, dt = dx/10

EDP Parabólica [ concepto ] [ ejercicio ]
Método explícito: [ Analítico ] [ Algoritmo ]
..

---

3. Algoritmo en Python

Se muestra una tabla resumida de resultados a forma de ejemplo.

Método explícito EDP Parabólica
lambda:  0.25
x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ] ... 1.0
Tabla de resultados en malla EDP Parabólica
j, U[:,: 10 ], primeras iteraciones
10 [60.   48.23 38.42 31.65 27.85 26.33 26.43 27.89 30.76 34.96 40.  ]
9 [60.   47.67 37.58 30.86 27.29 25.96 26.11 27.53 30.39 34.71 40.  ]
8 [60.   47.02 36.63 30.03 26.75 25.64 25.82 27.16 29.99 34.44 40.  ]
7 [60.   46.25 35.56 29.15 26.25 25.37 25.56 26.78 29.53 34.11 40.  ]
6 [60.   45.34 34.34 28.23 25.79 25.17 25.35 26.38 29.   33.72 40.  ]
5 [60.   44.21 32.93 27.29 25.41 25.05 25.18 25.98 28.4  33.23 40.  ]
4 [60.   42.77 31.29 26.37 25.14 25.   25.06 25.59 27.7  32.62 40.  ]
3 [60.   40.86 29.38 25.55 25.   25.   25.   25.23 26.88 31.8  40.  ]
2 [60.   38.12 27.19 25.   25.   25.   25.   25.   25.94 30.62 40.  ]
1 [60.   33.75 25.   25.   25.   25.   25.   25.   25.   28.75 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
>>> 

Se presenta la propuesta del algoritmo para el método explícito.

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
# Referencia: Chapra 30.2 p.888 pdf.912, Rodriguez 10.2 p.406
import numpy as np

# INGRESO
# Valores de frontera
Ta = 60
Tb = 40
#T0 = 25 # estado inicial de barra
fx = lambda x: 25 + 0*x # función inicial para T0
a = 0  # longitud en x
b = 1

K = 4     # Constante K
dx = 0.1  # Tamaño de paso
dt = dx/10

n = 100 # iteraciones en tiempo
verdigitos = 2   # decimales a mostrar en tabla de resultados

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
m = len(xi)
ultimox = m-1
ultimot = n-1
# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n+1), dtype=float)

# valores iniciales de u[:,j]
j = 0
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
tj = np.arange(0,(n+1)*dt,dt)
j = 0
while not(j>ultimot): # igual con lazo for
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    j=j+1

# SALIDA
np.set_printoptions(precision=verdigitos)
print('Método explícito EDP Parabólica')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:len(xi)],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',ultimox,'], primeras iteraciones')
for j in range(ultimox,-1,-1):
    print(j,u[:,j])

Si la cantidad de puntos aumenta al disminuir Δx y Δt, es mejor disminuir la cantidad de curvas a usar en el gráfico para evitar superponerlas. Se usa el parámetro ‘salto’ para indicar cada cuántas curvas calculadas se incorporan en la gráfica.

# GRAFICA ------------
import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneas
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
plt.show()

Note que en la gráfica se toma como coordenadas x vs t, mientras que para la solución de la malla en la matriz las se usan filas y columnas. En la matriz el primer índice es fila y el segundo índice es columna.

En el caso de una gráfica que se realice usando los ejes x,t,U en tres dimensiones se requiere usar las instrucciones:

# GRAFICA en 3D ------
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)

# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(m,tramos),dtype=float)
for k in range(0,tramos,1):
    U[:,k] = u[:,k*salto]
    tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)
U = np.transpose(U)

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[0],tk[1],U[1,0],'o',color ='orange')
graf_3D.plot(xi[1],tk[1],U[1,1],'o',color ='green')
graf_3D.plot(xi[2],tk[1],U[1,2],'o',color ='green')
graf_3D.plot(xi[1],tk[2],U[2,1],'o',color ='salmon',
             label='U[i,j+1]')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

EDP Parabólica [ concepto ] [ ejercicio ]
Método explícito: [ Analítico ] [ Algoritmo ]

EDP Parabólica [ concepto ] Método [ explícito ] [ implícito ]

..

---

1. EDP Parabólicas

Referencia:  Chapra 30.2 p.888 pdf.912, Burden 9Ed p714, Rodríguez 10.2 p406

Las Ecuaciones Diferenciales Parciales tipo parabólicas, semejantes a la mostrada, representa la ecuación de calor para una barra aislada sometida a fuentes de calor en cada extremo.

La temperatura se representa en el ejercicio como u[x,t]

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

o con vista en 3D:

Para la solución numérica, cambia la ecuación a su forma discreta usando diferencias finitas divididas que se sustituyen en la ecuación,

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

con lo que la ecuación continua se convierte a discreta:

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

Para interpretar mejor el resultado, se usa una malla que en cada nodo representa la temperatura como los valores u[x_i,t_j].

Para simplificar nomenclatura se usan los subíndices i para el eje de las x y j para el eje t, quedando u[ i , j ].

En el enunciado del problema habían establecido los valores en las fronteras:

• temperaturas en los extremos Ta, Tb
• la temperatura inicial de la barra T0,
• El parámetro para la barra K.

El resultado obtenido se interpreta como los valores de temperatura a lo largo de la barra luego de transcurrido un largo tiempo. Las temperaturas en los extremos de la barra varían entre Ta y Tb a lo largo del tiempo.

Tomando como referencia la malla, existirían algunas formas de plantear la solución, dependiendo de la diferencia finita usada: centrada, hacia adelante, hacia atrás.

EDP Parabólica [ concepto ] Método [ explícito ] [ implícito ]

---

But what is a partial differential equation? | DE2

3Blue1Brown. 2019 Abril 21

---

Tarea: Revisar ecuación para difusión de gases, segunda ley de Fick.

La difusión molecular desde un punto de vista microscópico y macroscópico.

Parabólicas: [ explícito ] [ implícito ]

Elípticas: [ iterativo ] [ implícito ]

Hiperbólicas