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} = 0La 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:
El término constante, por facilidad se simplifica con el valor de 1
\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} =1si 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 1para 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.