EDP y diferencias finitas

💡 En simple

Piensa en una plancha de metal con un borde caliente y los otros fríos. ¿Qué temperatura tiene cada punto cuando todo se estabiliza? Cada punto termina con el promedio de sus vecinos.

La computadora no puede trabajar con “infinitos puntos”, así que pone una cuadrícula (malla) y repite: “cada casilla = promedio de sus vecinas”, hasta que los números dejan de cambiar. Eso es resolver una EDP con diferencias finitas y el método de Jacobi.

Para el curso, lo importante es lo paralelo: cada casilla solo mira a sus vecinas. Si partes la malla entre procesos, cada uno solo necesita el borde del vecino (el halo de la U4.5).

Qué es una EDP

Una ecuación diferencial parcial (EDP) relaciona una función de varias variables, por ejemplo u(x, y) o u(x, t), con sus derivadas parciales. Modela calor, difusión, fluidos, ondas y potenciales gravitacionales o eléctricos (lámina 2).

Los tres tipos (láminas 5–6)

Cómo leerlo: la diferencia que importa en paralelo es si aparece el tiempo. Las elípticas no lo tienen: se itera hasta que la solución converge. Las parabólicas e hiperbólicas avanzan paso a paso en el tiempo, y ese avance no se puede paralelizar (lámina 31): solo se reparte el espacio.

Signo en la lámina 6

La lámina escribe la ecuación de onda como ∂²u/∂t² = −c² ∂²u/∂x², con signo menos. El signo correcto es +: ∂²u/∂t² = c² ∂²u/∂x². Con el signo menos la ecuación no describe una onda: es inestable y crece sin control. La lámina 35 ya la escribe con +.

De la derivada a la resta: diferencias finitas (láminas 4, 7–8)

Idea: la derivada es una pendiente. En una malla con separación h, la aproximas restando valores vecinos.

Derivada	Aproximación	Error
Primera (hacia adelante)	(ui+1 − ui) / h	O(h)
Primera (centrada)	(ui+1 − ui−1) / 2h	O(h²)
Segunda (centrada)	(ui+1 − 2ui + ui−1) / h²	O(h²)

Salen de la serie de Taylor: si sumas u(x₀+h) y u(x₀−h), los términos de orden impar se cancelan y queda u(x₀+h) + u(x₀−h) = 2u(x₀) + h²·u″(x₀) + O(h⁴). Al despejar u″ sale la fórmula de la segunda derivada.

Laplace 2D paso a paso (láminas 7–16)

Ecuación de Laplace (equilibrio de temperatura): ∂²u/∂x² + ∂²u/∂y² = 0.

1. Sustituir cada segunda derivada por su diferencia finita (con Δx = Δy = h):

   (u_i+1,j − 2u_i,j + u_i−1,j)/h² + (u_i,j+1 − 2u_i,j + u_i,j−1)/h² = 0

2. Multiplicar por h² y juntar términos:

   u_i+1,j + u_i−1,j + u_i,j+1 + u_i,j−1 − 4u_i,j = 0

3. Despejar u_i,j: u_i,j = ¼ · (suma de los 4 vecinos).

4. Frontera: si un vecino está en el borde, su valor es conocido y pasa al lado derecho (lámina 15): −4u_2,4 + u_3,4 + u_2,3 + u_2,5 = −u_1,4.

Una ecuación por punto interior: con N×N puntos hay N² ecuaciones con N² incógnitas (en 3D con 1000³ puntos, 10⁹ ecuaciones, lámina 24). Nadie resuelve eso por eliminación gaussiana: se itera.

Cómo leerlo: el punto naranja se recalcula con los 4 celestes. Este patrón fijo de vecinos se llama stencil de 5 puntos. Como cada punto solo depende de vecinos cercanos, el problema se reparte bien: un proceso solo necesita los puntos de sus bordes que tienen otros procesos.

Una observación sobre la lámina 9

La lámina escribe el paso como Δu = (promedio de vecinos) − uⁿ, y luego uⁿ⁺¹ = uⁿ + Δu. Las dos formas dan lo mismo: uⁿ⁺¹ = promedio de vecinos. El código de la lámina usa directamente la segunda.

Poisson: Laplace con fuente (láminas 23–25)

∇²u = S(x, y). La fuente S es la masa, la carga o una fuente de calor. Con el mismo procedimiento:

u^k+1_i,j = ¼ · ( u^k_i+1,j + u^k_i−1,j + u^k_i,j+1 + u^k_i,j−1 − h²·S_i,j )

Tipos de frontera: Dirichlet fija el valor (u = f en el borde). Neumann fija la derivada normal (∂u/∂n = g, por ejemplo “no sale calor”).

Jacobi vs Gauss-Seidel (láminas 25–26)

	Jacobi	Gauss-Seidel
Qué vecinos usa	todos de la iteración anterior (u_old)	los ya actualizados en esta iteración
Memoria	dos arreglos (u y u_old)	uno solo
Iteraciones para converger	más	menos (≈ la mitad)
Paralelismo	todos los puntos a la vez	el punto i espera al i−1: secuencial. Hay que reordenar (rojo-negro, conocimiento general)

Criterio de parada: diff = max |u^k+1 − u^k| < ε. En paralelo, cada proceso calcula su máximo local y luego se hace un MPI_Allreduce con MPI_MAX.

Complejidad: T(N) = k · m · N² (lámina 26)

• N: puntos por dimensión (N² en total en 2D).
• m: operaciones por punto (≈ 5 en el stencil: 3 sumas, 1 multiplicación y la diferencia del error).
• k: número de iteraciones hasta converger.

La trampa es que k no es constante: crece con N. Medido con laplace_sec.cpp (Laplace 1D, tolerancia 10⁻⁶):

N	Iteraciones de Jacobi hasta converger
50	6 066
100	21 980
200	77 614

Al duplicar N las iteraciones se multiplican por ≈ 3.6. Es decir, k crece casi como N². El motivo: el valor de la frontera avanza un punto por iteración, y además cada recorrido solo reduce un poco el error.

El código de clase no converge

laplace_sec.cpp usa max_iter = 10000, pero con N = 100 se necesitan 21 980 iteraciones. Al ejecutarlo se detiene en 10 000 con error 4.2·10⁻⁴, mil veces más que la tolerancia 10⁻⁶, y da u(0.4949) = 49.08 cuando el valor exacto es 49.49. La gráfica de clase “se ve” bien, pero no convergió. Si un examen pregunta por qué, la respuesta es que Jacobi converge lento: k ∝ N².

Práctica

1. Discretiza la ecuación del calor 1D ∂u/∂t = D·∂²u/∂x² con diferencia hacia adelante en el tiempo (paso Δt) y centrada en el espacio (paso h). Escribe u_iⁿ⁺¹.

Ver solución

Izquierda: (u_iⁿ⁺¹ − u_iⁿ)/Δt. Derecha: D·(u_i+1ⁿ − 2u_iⁿ + u_i−1ⁿ)/h².

Despejando: u_iⁿ⁺¹ = u_iⁿ + (D·Δt/h²)·(u_i+1ⁿ − 2u_iⁿ + u_i−1ⁿ).

Es el mismo stencil de 3 puntos, así que se paraleliza igual que Laplace 1D: 1 valor de halo por lado y por paso. Conocimiento general: este esquema explícito solo es estable si D·Δt/h² ≤ ½.

2. Una malla 3D de 512³ puntos, m = 8 operaciones por punto y k = 2000 iteraciones. ¿Cuántas operaciones hace Jacobi? Si la máquina hace 10 GFLOP/s, ¿cuánto tarda?

Ver solución

T = k·m·N³ (en 3D hay N³ puntos) = 2000 · 8 · 512³ = 2000 · 8 · 134 217 728 ≈ 2.15·10¹² operaciones.

A 10·10⁹ op/s: 2.15·10¹² / 10¹⁰ ≈ 215 s. Por eso estos problemas se reparten entre muchos procesos.

Autoevaluación

0/6 correctas · 0 respondidas Reiniciar

1. ¿Qué tipo de EDP es la ecuación de Laplace?

a) hiperbólicab) parabólicac) elípticad) ninguna

2. En una EDP parabólica u(x,t), ¿qué se puede paralelizar?

a) el tiempo y el espaciob) solo el tiempoc) solo el espacio (el tiempo es secuencial)d) nada

3. El stencil de Laplace 2D calcula u[i][j] como…

a) la suma de 8 vecinosb) ¼ de la suma de los 4 vecinosc) el máximo de los vecinosd) u_old[i][j] + 1

4. ¿Por qué Jacobi se paraleliza mejor que Gauss-Seidel?

a) porque converge más rápidob) porque solo usa valores de la iteración anterior, así que todos los puntos son independientesc) porque usa menos memoriad) porque no necesita frontera

5. En T(N) = k·m·N², ¿qué pasa con k al aumentar N en Jacobi?

a) queda constanteb) crece (≈ N²)c) disminuyed) depende solo de m

6. Una frontera de Neumann fija…

a) el valor de u en el bordeb) la derivada normal de u en el bordec) el número de iteracionesd) la fuente S

Fuente ·U5.1 Algoritmos-EDP.pdf (láminas 2–31) · U5.1 Algoritmos-EDP-code/laplace_sec.cpp