In [1]:
import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt
from ipywidgets import widgets
%matplotlib inline
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*}
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)
In [8]:
dp = 1
widgets.interact(plot,us=widgets.fixed(us1),i=(0,M,dp))
Out[8]:
In [9]:
us2 = Dirichlet(x,f1,g2,l,r,c**2,dx,dt,N,M)
In [10]:
dp = 1
widgets.interact(plot,us=widgets.fixed(us2),i=(0,M,dp))
Out[10]:
In [11]:
us3 = Dirichlet(x,f2,g2,l,r,c**2,dx,dt,N,M)
In [12]:
dp = 1
widgets.interact(plot,us=widgets.fixed(us3),i=(0,M,dp))
Out[12]:
In [ ]: