ILI286 - Computación Científica II

EDP Hiperbólicas: Diferencias Finitas

[S]cientific [C]omputing [T]eam

Version: 1.12

Librerías necesarias


In [1]:
import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt
from ipywidgets import widgets
%matplotlib inline

Ecuación de Onda 1D

Consideremos la siguiente EDP hiperbólica, para $ x \in [a,b] \ \ $ y $t \in [0,T_{max}] \ \ $: \begin{align*} u_{tt}(x,t) &= c^2\, u_{xx}(x,t) \\ u(x,0) &= f(x)\\ u_t(x,0) &= g(x)\\ u(a,t) &= l(t)\\ u(b,t) &= r(t) \end{align*}

Utilizando diferencias finitas, la ecuación nos queda: \begin{align*} u_{tt}(x,t) &= c^2\, u_{xx}(x,t) \\ \dfrac{u_i^{t+1}-2u_i^t+u_i^{t-1}}{\Delta t^2} &= c^2\, \dfrac{u_{i+1}^t-2u_i^t+u_{i-1}^t}{\Delta x^2} \end{align*}

Nuestra incógnita sería $u_i^{t+1}$, la cual se puede calcular teniendo información de los dos tiempos anteriores. Por lo tanto, despejando llegamos a:

$$u_i^{t+1} = \dfrac{c^2 \Delta t^2}{\Delta x^2}\left(u_{i+1}^t + u_{i-1}^t\right) + 2\left(1 - \dfrac{c^2 \Delta t^2}{\Delta x^2}\right)u_i^t - u_i^{t-1}$$

Si bien esto se resuelve directamente a medida que avanzamos en el tiempo, en $t=1$ nos faltará el valor de $u_i^{-1}$, pero recordando que: $$u_t = \dfrac{u_i^{t+1} - u_i^{t-1}}{2\Delta t}$$ se puede utilizar la condición incial para obtener: \begin{align*} u_i^{t+1} &= \dfrac{c^2 \Delta t^2}{2 \Delta x^2}\left(u_{i+1}^t + u_{i-1}^t\right) + \left(1 - \dfrac{c^2 \Delta t^2}{\Delta x^2}\right)u_i^t + \Delta t u_t \\ u_i^1 &= \dfrac{c^2 \Delta t^2}{2 \Delta x^2}\left(u_{i+1}^0 + u_{i-1}^0\right) + \left(1 - \dfrac{c^2 \Delta t^2}{\Delta x^2}\right)u_i^0 + \Delta t u_t(x_i,0) \end{align*}

Finite Differences con condición de borde de Dirichlet


In [2]:
def Dirichlet(x,f,g,l,r,c2,dx,dt,N,M):
    
    sigma2 = c2 * dt**2 / dx**2
    
    print('sigma**2= ',sigma2)
        
    # The approximation for all time-steps
    w  = np.zeros((M+1,N))
    
    # Initial condition
    w[0,:]=f
        
    # Building A
    d=np.zeros(N-2)
    d[0]=2-2*sigma2
    d[1]=sigma2
    A=spla.toeplitz(d)
    
    ## First step
    tmp=np.zeros(N-2)
    tmp[0]=w[0,0]
    tmp[-1]=w[0,-1]
    w[1,1:-1]=0.5*(np.dot(A,w[0,1:-1]))+dt*g[1:-1]+0.5*sigma2*tmp
    
    ## All the rest of the steps
    for t in range(2,M+1):
        tmp[0]=l(dt*t)
        tmp[-1]=r(dt*t)
        w[t,1:-1]=np.dot(A,w[t-1,1:-1])-w[t-2,1:-1]+sigma2*tmp
        w[t,0]=l(dt*t)
        w[t,-1]=r(dt*t)
    return w

In [3]:
def l(t):
    return 0*t

def r(t):
    return 0*t

In [4]:
def plot(us,i=0):
    plt.plot(x,us[i])
    plt.ylim(-0.11,0.11)
    plt.title('Tiempo: '+str(i*dt))
    plt.grid()
    plt.show()

In [5]:
N=100
M=200
TMAX = 40
x = np.linspace(-20,20,N)
t = np.linspace(0,TMAX,M)

c = 1

dx = x[1]-x[0]
dt = t[1]-t[0]

# RECALL: CFL condition:
#  if c dt/dx<=1 => stability
#  else (c dt/dx>1) instability.

In [6]:
## Numerical Initial Conditions, notice these are numerical arrays!
# i.e. vectors not functions.
# u(x,0)=f(x)
f1 = 0.1*np.exp(-x**2)
f2 = 0.1*np.sin(2*x*np.pi/8)

# u_t(x,0)=f(x)
# SOLITON
g1 = (np.exp(-(x-0.2)**2) - np.exp(-(x+0.2)**2))/2
# ONDA QUE SE DIVIDE
g2 = np.zeros(N)

In [7]:
us1 = Dirichlet(x,f1,g1,l,r,c**2,dx,dt,N,M)


sigma**2=  0.24749375015782332

In [8]:
dp = 1
widgets.interact(plot,us=widgets.fixed(us1),i=(0,M,dp))


Out[8]:
<function __main__.plot(us, i=0)>

In [9]:
us2 = Dirichlet(x,f1,g2,l,r,c**2,dx,dt,N,M)


sigma**2=  0.24749375015782332

In [10]:
dp = 1
widgets.interact(plot,us=widgets.fixed(us2),i=(0,M,dp))


Out[10]:
<function __main__.plot(us, i=0)>

In [11]:
us3 = Dirichlet(x,f2,g2,l,r,c**2,dx,dt,N,M)


sigma**2=  0.24749375015782332

In [12]:
dp = 1
widgets.interact(plot,us=widgets.fixed(us3),i=(0,M,dp))


Out[12]:
<function __main__.plot(us, i=0)>

Acknowledgements

  • Material creado por profesor Claudio Torres (ctorres@inf.utfsm.cl) y ayudantes: Alvaro Salinas y Martín Villanueva. DI UTFSM. Noviembre 2016.

[Update 2019] (C. Torres) Fixing titles and small changes.


In [ ]: