Ejercicios resueltos de diferencias finitas

Dividimos $[0{,}1] \times [0{,}1]$ en una malla $3 \times 3$:

$h = \frac{1}{3}$

Nodos: $(x_i, y_j)$ donde $x_i = ih$ y $y_j = jh$ para $i,j = 0{,}1{,}2{,}3$

Los nodos interiores son: $(1/3, 1/3)$, $(1/3, 2/3)$, $(2/3, 1/3)$, $(2/3, 2/3)$

$\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2}$

$\frac{\partial^2 u}{\partial y^2} \approx \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2}$

La ecuación de Poisson se convierte en:

$\frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2} = -x_i y_j$

$u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j} = -h^2 x_i y_j$

Sean $u_1 = u_{1{,}1}$, $u_2 = u_{1{,}2}$, $u_3 = u_{2{,}1}$, $u_4 = u_{2{,}2}$

Nodo $(1/3, 1/3)$

$u_{2{,}1} + u_{0{,}1} + u_{1{,}2} + u_{1{,}0} - 4u_{1{,}1} = -h^2 \cdot \frac{1}{3} \cdot \frac{1}{3}$

Como $u_{0{,}1} = u_{1{,}0} = 0$ (condiciones de frontera):

$u_3 + 0 + u_2 + 0 - 4u_1 = -\frac{1}{81}$

$u_3 + u_2 - 4u_1 = -\frac{1}{81}$

Nodo $(1/3, 2/3)$

$u_4 + 0 + 0 + u_1 - 4u_2 = -h^2 \cdot \frac{1}{3} \cdot \frac{2}{3} = -\frac{2}{81}$

$u_4 + u_1 - 4u_2 = -\frac{2}{81}$

Nodo $(2/3, 1/3)$

$0 + u_1 + u_4 + 0 - 4u_3 = -h^2 \cdot \frac{2}{3} \cdot \frac{1}{3} = -\frac{2}{81}$

$u_1 + u_4 - 4u_3 = -\frac{2}{81}$

Nodo $(2/3, 2/3)$

$0 + u_2 + 0 + u_3 - 4u_4 = -h^2 \cdot \frac{2}{3} \cdot \frac{2}{3} = -\frac{4}{81}$

$u_2 + u_3 - 4u_4 = -\frac{4}{81}$

$\begin{bmatrix} -4 & 1 & 1 & 0 \\ 1 & -4 & 0 & 1 \\ 1 & 0 & -4 & 1 \\ 0 & 1 & 1 & -4 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} = \begin{bmatrix} -1/81 \\ -2/81 \\ -2/81 \\ -4/81 \end{bmatrix}$

Usando eliminación gaussiana o métodos iterativos:

$u_1 \approx 0.00977$

$u_2 \approx 0.01337$

$u_3 \approx 0.01337$

$u_4 \approx 0.01903$

Resultado:

$u(1/3, 1/3) \approx 0.00977$

$u(1/3, 2/3) \approx 0.01337$

$u(2/3, 1/3) \approx 0.01337$

$u(2/3, 2/3) \approx 0.01903$

Espacial: $x_i = ih$ donde $h = \frac{1}{N}$ y $i = 0, 1, 2, ..., N$

Temporal: $t_n = n\Delta t$ donde $n = 0, 1, 2, ...$

Sea $u_i^n = u(x_i, t_n)$ la aproximación en el nodo $(x_i, t_n)$.

$\frac{\partial u}{\partial t} \approx \frac{u_i^{n+1} - u_i^n}{\Delta t}$

$\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{h^2}$

$\frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{h^2}$

Definiendo $r = \frac{\Delta t}{h^2}$:

$u_i^{n+1} - u_i^n = r(u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1})$

$-ru_{i-1}^{n+1} + (1 + 2r)u_i^{n+1} - ru_{i+1}^{n+1} = u_i^n$

Para $i = 1, 2, ..., N-1$ (nodos interiores):

$\begin{bmatrix} 1+2r & -r & 0 & 0 & \cdots & 0 \\ -r & 1+2r & -r & 0 & \cdots & 0 \\ 0 & -r & 1+2r & -r & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & -r & 1+2r \end{bmatrix} \begin{bmatrix} u_1^{n+1} \\ u_2^{n+1} \\ u_3^{n+1} \\ \vdots \\ u_{N-1}^{n+1} \end{bmatrix} = \begin{bmatrix} u_1^n \\ u_2^n \\ u_3^n \\ \vdots \\ u_{N-1}^n \end{bmatrix}$

Condiciones de frontera: $u_0^{n+1} = u_N^{n+1} = 0$ para todo $n$

Condición inicial: $u_i^0 = \sin(\pi x_i) = \sin(\pi ih)$ para $i = 0, 1, ..., N$

Para n = 0, 1, 2, ... (pasos temporales):

1. Formar vector RHS = $[u_1^n, u_2^n, \ldots, u_{N-1}^n]^T$
2. Resolver sistema tridiagonal $Au^{n+1} = u^n$
3. Aplicar condiciones de frontera: $u_0^{n+1} = u_N^{n+1} = 0$

$h = 1/4 = 0.25$, $r = \frac{0.01}{(0.25)^2} = 0.16$

Matriz del sistema:

$\begin{bmatrix} 1.32 & -0.16 & 0 \\ -0.16 & 1.32 & -0.16 \\ 0 & -0.16 & 1.32 \end{bmatrix}$

Condición inicial:

$u_1^0 = \sin(\pi/4) = \sqrt{2}/2 \approx 0.707$

$u_2^0 = \sin(\pi/2) = 1$

$u_3^0 = \sin(3\pi/4) = \sqrt{2}/2 \approx 0.707$

1. Estabilidad incondicional: No hay restricción en $\Delta t$ para estabilidad
2. Precisión: Error de truncamiento $O(\Delta t + h^2)$
3. Robustez: Permite pasos temporales grandes

Desventaja: Requiere resolver sistema lineal en cada paso temporal.

Para $\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$, el esquema de Lax-Wendroff es:

$u_j^{n+1} = u_j^n - \frac{c\Delta t}{2h}(u_{j+1}^n - u_{j-1}^n) + \frac{c^2(\Delta t)^2}{2h^2}(u_{j+1}^n - 2u_j^n + u_{j-1}^n)$

donde $h = \Delta x$ es el espaciado espacial.

Sea $\nu = \frac{c\Delta t}{h}$ el número de Courant.

El esquema se convierte en:

$u_j^{n+1} = u_j^n - \frac{\nu}{2}(u_{j+1}^n - u_{j-1}^n) + \frac{\nu^2}{2}(u_{j+1}^n - 2u_j^n + u_{j-1}^n)$

Supongamos una solución de la forma $u_j^n = \xi^n e^{ikjh}$ donde:

• $\xi$ es el factor de amplificación
• $k$ es el número de onda
• $i = \sqrt{-1}$

$\xi^{n+1} e^{ikjh} = \xi^n e^{ikjh} - \frac{\nu}{2}\xi^n(e^{ik(j+1)h} - e^{ik(j-1)h})$

$+ \frac{\nu^2}{2}\xi^n(e^{ik(j+1)h} - 2e^{ikjh} + e^{ik(j-1)h})$

$\xi = 1 - \frac{\nu}{2}(e^{ikh} - e^{-ikh}) + \frac{\nu^2}{2}(e^{ikh} - 2 + e^{-ikh})$

Usando $e^{ikh} - e^{-ikh} = 2i\sin(kh)$ y $e^{ikh} + e^{-ikh} = 2\cos(kh)$:

$\xi = 1 - \frac{\nu}{2}(2i\sin(kh)) + \frac{\nu^2}{2}(2\cos(kh) - 2)$

$\xi = 1 - i\nu\sin(kh) + \nu^2(\cos(kh) - 1)$

$\cos(kh) - 1 = -2\sin^2(kh/2)$

$\xi = 1 - i\nu\sin(kh) - 2\nu^2\sin^2(kh/2)$

$|\xi|^2 = \xi \bar{\xi} = (1 - 2\nu^2\sin^2(kh/2))^2 + \nu^2\sin^2(kh)$

Usando $\sin^2(kh) = 4\sin^2(kh/2)\cos^2(kh/2)$:

$|\xi|^2 = (1 - 2\nu^2\sin^2(kh/2))^2 + 4\nu^2\sin^2(kh/2)\cos^2(kh/2)$

$|\xi|^2 = 1 - 4\nu^2\sin^2(kh/2) + 4\nu^4\sin^4(kh/2) + 4\nu^2\sin^2(kh/2)\cos^2(kh/2)$

$= 1 - 4\nu^2\sin^2(kh/2) + 4\nu^2\sin^2(kh/2)[\nu^2\sin^2(kh/2) + \cos^2(kh/2)]$

$= 1 - 4\nu^2\sin^2(kh/2)(1 - \cos^2(kh/2) - \nu^2\sin^2(kh/2))$

$= 1 - 4\nu^2\sin^2(kh/2)\sin^2(kh/2)(1 - \nu^2)$

$= 1 - 4\nu^2\sin^4(kh/2)(1 - \nu^2)$

Para estabilidad, necesitamos $|\xi|^2 \leq 1$ para todo $k$.

$1 - 4\nu^2\sin^4(kh/2)(1 - \nu^2) \leq 1$

Esto requiere $4\nu^2\sin^4(kh/2)(1 - \nu^2) \geq 0$.

Como $\sin^4(kh/2) \geq 0$ y $\nu^2 \geq 0$, necesitamos $(1 - \nu^2) \geq 0$.

Condición de estabilidad: $\nu^2 \leq 1$, o equivalentemente:

$$\frac{c\Delta t}{h} \leq 1$$

Esta es la condición CFL (Courant-Friedrichs-Lewy) para el esquema de Lax-Wendroff.

Interpretación física: El paso temporal debe ser lo suficientemente pequeño para que la información no viaje más de una celda espacial por paso temporal.

Para obtener una aproximación de cuarto orden de $\frac{\partial u}{\partial x}$, usamos 5 puntos:

$\frac{\partial u}{\partial x}\Big|_{x_j} \approx \frac{-u_{j+2} + 8u_{j+1} - 8u_{j-1} + u_{j-2}}{12h} + O(h^4)$

$u_{j+1} = u_j + h\frac{\partial u}{\partial x} + \frac{h^2}{2!}\frac{\partial^2 u}{\partial x^2} + \frac{h^3}{3!}\frac{\partial^3 u}{\partial x^3} + \frac{h^4}{4!}\frac{\partial^4 u}{\partial x^4} + \frac{h^5}{5!}\frac{\partial^5 u}{\partial x^5} + O(h^6)$

$u_{j-1} = u_j - h\frac{\partial u}{\partial x} + \frac{h^2}{2!}\frac{\partial^2 u}{\partial x^2} - \frac{h^3}{3!}\frac{\partial^3 u}{\partial x^3} + \frac{h^4}{4!}\frac{\partial^4 u}{\partial x^4} - \frac{h^5}{5!}\frac{\partial^5 u}{\partial x^5} + O(h^6)$

$u_{j+2} = u_j + 2h\frac{\partial u}{\partial x} + \frac{4h^2}{2!}\frac{\partial^2 u}{\partial x^2} + \frac{8h^3}{3!}\frac{\partial^3 u}{\partial x^3} + \frac{16h^4}{4!}\frac{\partial^4 u}{\partial x^4} + \frac{32h^5}{5!}\frac{\partial^5 u}{\partial x^5} + O(h^6)$

$u_{j-2} = u_j - 2h\frac{\partial u}{\partial x} + \frac{4h^2}{2!}\frac{\partial^2 u}{\partial x^2} - \frac{8h^3}{3!}\frac{\partial^3 u}{\partial x^3} + \frac{16h^4}{4!}\frac{\partial^4 u}{\partial x^4} - \frac{32h^5}{5!}\frac{\partial^5 u}{\partial x^5} + O(h^6)$

Buscamos coeficientes $a, b, c, d$ tales que:

$a u_{j-2} + b u_{j-1} + c u_{j+1} + d u_{j+2} = A h \frac{\partial u}{\partial x} + O(h^5)$

donde queremos $A = 1$ y que se cancelen los términos de $h^2$, $h^3$, y $h^4$.

Condiciones:

• Coeficiente de $u_j$: $a + b + c + d = 0$
• Coeficiente de $h\frac{\partial u}{\partial x}$: $-2a - b + c + 2d = A = 1$
• Coeficiente de $h^2\frac{\partial^2 u}{\partial x^2}$: $2a + \frac{b}{2} + \frac{c}{2} + 2d = 0$
• Coeficiente de $h^3\frac{\partial^3 u}{\partial x^3}$: $-\frac{4a}{3} - \frac{b}{6} + \frac{c}{6} + \frac{4d}{3} = 0$

Resolviendo: $a = \frac{1}{12}$, $b = -\frac{8}{12}$, $c = \frac{8}{12}$, $d = -\frac{1}{12}$

Por tanto:

$\frac{\partial u}{\partial x} \approx \frac{u_{j-2} - 8u_{j-1} + 8u_{j+1} - u_{j+2}}{12h}$

$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$

Usando diferencias hacia adelante en tiempo (Euler explícito):

$\frac{u_j^{n+1} - u_j^n}{\Delta t} + c \frac{u_{j-2}^n - 8u_{j-1}^n + 8u_{j+1}^n - u_{j+2}^n}{12h} = 0$

$u_j^{n+1} = u_j^n - \frac{c\Delta t}{12h}(u_{j-2}^n - 8u_{j-1}^n + 8u_{j+1}^n - u_{j+2}^n)$

Substituyendo $u_j^n = \xi^n e^{ikjh}$:

$\xi = 1 - \frac{c\Delta t}{12h}(e^{-2ikh} - 8e^{-ikh} + 8e^{ikh} - e^{2ikh})$

$= 1 - \frac{c\Delta t}{12h} \cdot 2i[8\sin(kh) - \sin(2kh)]$

$= 1 - i\frac{c\Delta t}{6h}[8\sin(kh) - \sin(2kh)]$

$|\xi|^2 = 1 + \left(\frac{c\Delta t}{6h}\right)^2[8\sin(kh) - \sin(2kh)]^2 \geq 1$

La desigualdad es estricta para casi todos los modos $k$, sea cual sea $\Delta t$: con Euler explícito en el tiempo el esquema es incondicionalmente inestable (igual que con la diferencia centrada de segundo orden).

Solución práctica: avanzar en el tiempo con Runge-Kutta de cuarto orden (RK4). Entonces el esquema es estable si $\frac{c\Delta t}{h} \leq 2.06$ aproximadamente (frente a $2.83$ con la diferencia centrada de segundo orden).

Para j = 2, 3, ..., N-3:

$u_j^{n+1} = u_j^n - (cdt)$/(12h)

$(u[j-2]^n - 8u[j-1]^n + 8u[j+1]^n - u[j+2]^n)$

Para los puntos cerca de la frontera (j = 0, 1, N-2, N-1):

Usar esquemas de menor orden o condiciones de frontera especiales

Ventajas:

• Mayor precisión: Error $O(h^4)$ vs $O(h^2)$ para esquemas estándar
• Menor dispersión numérica para ondas de alta frecuencia

Desventajas:

• Más restrictiva en estabilidad
• Requiere más puntos (complicaciones en fronteras)
• Mayor costo computacional por punto

Aplicación: Especialmente útil para problemas de propagación de ondas donde se requiere alta precisión en largas distancias.