La PizarraMatemática.

Matemáticas para ingeniería

Diferencias finitas y discretización de ecuaciones

Si no podemos resolver una ecuación para cada punto de un intervalo, podemos empezar calculando valores en una malla. La idea parece sencilla: sustituir derivadas por diferencias de vecinos. Lo haremos despacio, cuidando el paso, el signo y las fronteras. Después comprobaremos por qué una malla más fina ayuda en ciertos casos y puede fallar en otros.

1. La malla define los puntos y las incógnitas que de verdad quedan

Dividir [a, b] en N intervalos iguales da h=(b−a)/N y nodos xᵢ=a+ih, i=0,…, N.

Son N+1 nodos. Con valores de Dirichlet fijados en ambos extremos quedan N−1 incógnitas interiores. El vector U aproxima valores de una función; aún no define por sí solo qué pasa entre nodos.

Contar intervalos, nodos e incógnitas evita cambiar el sistema
[0,1]; N=4; h=1/4U00U11/4U21/2U33/4U414 intervalos; 5 nodos3 incógnitas interioresU₀ y U₄ vienen de la frontera.Una fila para cada nodo interior.
El intervalo cero a uno se divide en cuatro partes iguales, con cinco nodos. Los extremos U cero y U cuatro están fijados por condiciones de Dirichlet; quedan tres incógnitas U uno, U dos y U tres. El paso es un cuarto y cada incógnita interior aporta una ecuación.

Las fórmulas de paso uniforme necesitan distancias iguales. En una malla no uniforme hay que volver a usar sus anchuras reales. No se elige una h cómoda y se la aplica a vecinos que están a otras distancias.

Contar antes de construir la matriz

Intervalo [0,1], N=4 y fronteras conocidas.

  1. Paso

    h=(1−0)/4=1/4.
  2. Nodos

    Son 0, 1/4, 1/2, 3/4 y uno: cinco en total.

  3. Incógnitas

    U₀ y U₄ están fijadas. Se resuelve para U₁, U₂ y U₃: tres ecuaciones interiores.

Resultado

Cuatro intervalos no son cuatro incógnitas interiores.

Lo comprobamos

Si N=2 hay tres nodos y solo una incógnita interior.

Una anchura desigual cambia la fórmula

Malla 0,1,3; u(x)=x². Aproximar u″ en x=1.

  1. Distancias

    A la izquierda h₋=1; a la derecha h₊=2.
    Las pendientes de las cuerdas son (1−0)/1=1 y (9−1)/2=4.
  2. Usar ambas

    La fórmula 2[(u derecha−u centro)/h₊−(u centro−u izquierda)/h₋]/(h₋+h₊) da 2(4−1)/3=2.
  3. Error de malla

    Usar como si todo tuviera h=1 daría 0−2·1+9=7.

    La fórmula no uniforme es exacta para esta cuadrática; su orden general requiere analizar cómo se refinan ambas distancias.

Resultado

El valor dos coincide con u″, usando la geometría correcta.

Lo comprobamos

No se promete orden dos en cualquier familia arbitraria de mallas desiguales.

2. Las diferencias de primera derivada comparan pendientes y cancelan errores distintos

La diferencia hacia delante es [u(x+h)−u(x)]/h; hacia atrás, [u(x)−u(x−h)]/h. Si u tiene segunda derivada acotada, Taylor muestra error O(h). La centrada [u(x+h)−u(x−h)]/(2h) cancela el término de segundo grado y tiene error O(h²) con tercera derivada acotada.

La suavidad es una hipótesis real. Por ejemplo, la fórmula centrada de |x| en cero da cero para cualquier h, pero |x| no es derivable allí. Una operación numérica definida no demuestra que exista el objeto que pretendemos aproximar.

Tres pendientes para una parábola

u(x)=x², x=2 y h=1/2; u′(2)=4.

  1. Delante

    [(5/2)²−2²]/(1/2)=(25/4−4)/(1/2)=9/2.
  2. Detrás

    [2²−(3/2)²]/(1/2)=(4−9/4)/(1/2)=7/2.
  3. Centrada

    [(5/2)²−(3/2)²]/1=4.

    Las dos cuerdas laterales fallan por ±1/2; al centrar se cancela ese error para esta cuadrática.

Resultado

La centrada es exacta en este ejemplo especial.

Lo comprobamos

Esa exactitud no significa que todas las funciones tengan error cero.

Un cúbico deja ver el orden de las fórmulas

u(x)=x³ en x=1, con h>0; la derivada exacta es tres.

  1. Expandir delante

    [(1+h)³−1]/h=3+3h+h².

    El error es 3h+h².

  2. Expandir centrada

    [(1+h)³−(1−h)³]/(2h)=3+h².

    El error es h².

  3. Números

    Con h=0,1 dan 3,31 y 3,01.

    Con h=0,05, el error centrado pasa de 0,01 a 0,0025, una cuarta parte.

Resultado

El orden describe cómo baja el error, no cuántos decimales imprimimos.

Lo comprobamos

Las potencias se expanden antes de cancelar; el denominador centrado es 2h.

3. La segunda diferencia necesita dividir por el cuadrado del paso

Taylor a ambos lados da u(x±h)=u(x)±hu′(x)+(h²/2)u″(x)±(h³/6)u‴(x)+O(h⁴), con la regularidad pertinente.

Sumar y restar 2u(x) cancela términos impares: u(x−h)−2u(x)+u(x+h)=h²u″(x)+O(h⁴).

Dividir por h² produce una aproximación de orden dos si la cuarta derivada está controlada. Sin h² cambiamos escala y unidades. Los errores de los valores usados también se amplifican por esa división; no todo refinamiento beneficia a datos ruidosos.

Un cuarto grado deja un error que podemos calcular exactamente

u(x)=x⁴ en x=1, h=1/2; u″(1)=12.

  1. Valores

    u(1/2)=1/16, u(1)=1 y u(3/2)=81/16.
  2. Combinar

    1/16−2+81/16=25/8.
  3. Dividir

    (25/8)/(1/4)=25/2=12,5.

    Para esta función la fórmula general da 12+2h²; el error es 1/2.

Resultado

La fórmula aproxima curvatura, no solo suma tres alturas.

Lo comprobamos

Con h=1/4 el error sería 1/8, una cuarta parte del anterior.

Las perturbaciones pequeñas de altura no desaparecen al dividir

Cada uno de los tres valores tiene error de módulo como máximo ε.

  1. Numerador

    El error combinado es η izquierda−2η centro+η derecha; su módulo es como máximo 4ε.

  2. Curvatura

    Al dividir por h², la perturbación queda acotada por 4ε/h².

  3. Refinar

    Reducir h a la mitad multiplica esa cota por cuatro. El error de truncamiento baja, pero esta otra fuente puede subir.

Resultado

Hay que controlar suavidad y calidad de los datos.

Lo comprobamos

La cota es una garantía de peor caso; las perturbaciones concretas pueden cancelar parte de su efecto.

4. Cada nodo interior aporta una ecuación con el signo del operador

Para −u″=f, cambiamos el signo de la segunda diferencia: (2Uᵢ−Uᵢ₋₁−Uᵢ₊₁)/h²=f(xᵢ).

La diagonal es positiva y los vecinos negativos. Al multiplicar por h² también multiplicamos el segundo miembro: 2Uᵢ−Uᵢ₋₁−Uᵢ₊₁=h²f(xᵢ).

Las dos formas son equivalentes si se escala todo. Mezclarlas altera el problema. En una dimensión, cada incógnita solo se conecta a sus dos vecinas y aparece una matriz tridiagonal.

Dos nodos interiores se convierten en dos filas

−u″=6 en [0,1], u(0)=u(1)=0, N=3.

  1. Paso

    h=1/3; interiores en 1/3 y 2/3.
  2. Sistema

    Sin dividir por h²: 2U₁−U₂=2/3 y −U₁+2U₂=2/3.

    En la otra escala la matriz tiene diagonal 18 y vecinos −9, con segundo miembro seis.

  3. Resolver

    Por simetría U₁=U₂, y cada uno vale 2/3.

    La solución exacta 3x(1−x) da esos mismos valores.

Resultado

El signo de menos u″ produce una matriz con diagonal positiva.

Lo comprobamos

La coincidencia nodal aquí se relaciona con la exactitud de la segunda diferencia para cuadráticas.

Una solución no cuadrática ya muestra error de discretización

−u″=12x², fronteras cero en [0,1], h=1/2.

  1. Nodo

    Solo x=1/2 es interior, y f allí vale tres.
  2. Fila

    2U/(1/2)²=3, así 8U=3 y U=3/8.
  3. Comparar

    La función exacta u=x−x⁴ satisface los datos y vale 7/16 en el centro. El error numérico U−u es −1/16, aunque la fila 8U=3 se cumpla exactamente.

Resultado

Resolver el sistema discreto exactamente no elimina el error de discretización.

Lo comprobamos

El residual de esa fila y el error respecto de la función continua son cantidades diferentes.

5. Las fronteras conocidas pasan a las filas próximas y Neumann plantea otra tarea

Con U₀=α y U_N=β, sus términos se trasladan al segundo miembro de las filas vecinas. Para la forma sin dividir por h², se añade α a la primera carga y β a la última. No se vuelven a tratar como incógnitas sin añadir sus condiciones.

Neumann fija una derivada en lugar de un valor y requiere otra discretización del borde. En un problema puramente Neumann, una constante puede quedar indeterminada y el dato f necesita compatibilidad. No se conserva sin cambios la matriz de Dirichlet.

Una frontera no nula aparece aunque f sea cero

−u″=0, u(0)=1, u(1)=3, N=4.

  1. Filas

    Para U₁, U₂, U₃, la matriz tiene diagonal dos y vecinos menos uno, y la carga es (1,0,3).

  2. Solución

    U=(3/2;2;5/2).
    Por ejemplo, 2·(3/2)−2=1 y −2+2·(5/2)=3.
  3. Interpretar

    Son valores de la recta 1+2x. Un segundo miembro nulo habría perdido las dos condiciones de frontera.

Resultado

El dato del borde también es parte de la ecuación cercana.

Lo comprobamos

En la escala dividida por h², la carga de las fronteras se divide igualmente por h².

Una condición de compatibilidad puede impedir que exista solución

−u″=f en [0,1], con u′(0)=u′(1)=0.

  1. Integrar

    ∫₀¹f dx=u′(0)−u′(1)=0.
  2. Dato incompatible

    Si f=1, la integral vale uno.

    No existe una solución que satisfaga ambas derivadas.

  3. Constante libre

    Si f=0, cualquier constante satisface la ecuación y ambas condiciones. Fijar además media cero selecciona la constante cero.

Resultado

Existencia y unicidad cambian con el tipo de frontera.

Lo comprobamos

No se fuerza una respuesta única eliminando una fila sin justificar qué problema queda.

6. Consistencia y estabilidad conectan el defecto local con el error global

El defecto de truncamiento se obtiene sustituyendo valores de la función exacta en el operador discreto Lₕ. El error global compara la solución U con esos valores exactos.

Si τ=Lₕu−f, entonces Lₕ(U−u)=−τ.

Para deducir error pequeño necesitamos controlar Lₕ⁻¹ con la escala y norma correctas.

Para este Dirichlet unidimensional en [0,1], el principio de comparación de la tridiagonal permite la cota ‖Lₕ⁻¹g‖∞≤‖g‖∞/8.

Una barrera es w=x(1−x)/2: Lₕw=1 y 0≤w≤1/8.

Las barreras ±‖g‖∞w acotan la solución. Esta constante corresponde al operador que incluye 1/h², no a cualquier matriz reescalada.

Una cota estable convierte truncamiento de orden dos en error de orden dos

u tiene cuarta derivada acotada por M y satisface fronteras de Dirichlet.

  1. Consistencia

    Taylor para la fórmula centrada da ‖τ‖∞≤Mh²/12 con las hipótesis suaves habituales.

  2. Estabilidad

    Aplicar la cota del inverso da ‖U−u‖∞≤(1/8)(Mh²/12)=Mh²/96.
  3. Significado

    La conclusión combina ambas propiedades. Una fórmula consistente sin control del inverso no daría por sí sola esta garantía.

Resultado

La precisión global se justifica con condiciones explícitas.

Lo comprobamos

En la comparación se usa la misma función en los mismos nodos y fronteras exactas.

Un ejemplo permite verificar la cota sin conocer U por tanteo

u=x−x⁴, f=12x² y extremos cero, con malla uniforme.

  1. Defecto

    Como la segunda diferencia de x⁴ añade 2h², Lₕu=f+2h².
    Así τ=2h².
  2. Error exacto

    La función nodal e=−h²x(1−x) tiene fronteras cero y Lₕe=−2h².
    Por unicidad discreta, e=U−u.
  3. Máximo

    Si N es par, el centro es un nodo y ‖e‖∞=h²/4.

    Coincide con ‖τ‖∞/8; para N impar el máximo nodal es menor.

Resultado

Para N par, h=1/2,1/4,1/8 produce errores 1/16,1/64,1/256.

Lo comprobamos

La cota de comparación explica a la vez el orden y la constante de este caso.

7. En difusión, refinar espacio obliga a revisar también el paso de tiempo

Para uₜ=D uₓₓ, D>0, Euler explícito en tiempo y diferencia centrada en espacio dan Uᵢ nuevo=λUᵢ₋₁+(1−2λ)Uᵢ+λUᵢ₊₁, con λ=DΔt/h². Si 0≤λ≤1/2, los pesos no son negativos y suman uno. Esto da una garantía de estabilidad en el marco unidimensional usual con fronteras controladas.

El paso estable combina vecinos con pesos no negativos
Difusión explícita: λ=2/5Pesos sobre el nivel anterior2402/51/52/5U nuevo = 8/5Entre el mínimo 0 y el máximo 4.Pesos no negativos; suma uno.Con λ>1/2, esta garantía se pierde.
Para lambda dos quintos, los pesos de izquierda, centro y derecha son dos quintos, un quinto y dos quintos. Con valores dos, cuatro y cero, la salida es ocho quintos, dentro del rango anterior. La condición lambda no mayor que un medio asegura esa combinación convexa en esta difusión explícita unidimensional.

Para una onda discreta, el factor de amplificación es 1−4λsin²(θ/2). Fuera del intervalo permitido algunas perturbaciones pueden crecer. Euler implícito requiere resolver un sistema por paso y tiene otro factor: 1/[1+4λsin²(θ/2)]. Su estabilidad para λ≥0 en este problema lineal no elimina el error temporal.

Ver los pesos de un paso y una perturbación que puede crecer

Tomar λ=2/5 y valores antiguos izquierdo dos, central cuatro, derecho cero.

  1. Combinar

    U nuevo=(2/5)·2+(1/5)·4+(2/5)·0=8/5, dentro del rango cero a cuatro.

  2. Cambiar λ

    En una malla infinita o periódica con alternancia de signos, θ=π.

    Con λ=3/5, su factor es 1−4λ=−7/5.
  3. Efecto

    Cada paso invierte el signo y multiplica la amplitud por 1,4. Esa perturbación crece, aunque la difusión continua deba amortiguarla.

Resultado

La restricción de paso evita una amplificación numérica incompatible con el modelo.

Lo comprobamos

El ejemplo de alternancia declara su malla; no se impone sin revisión a cualquier frontera finita.

La misma Δt puede dejar de ser segura al afinar h

D=1, h=1/10 y Δt=1/250.

  1. Antes

    λ=(1/250)/(1/100)=2/5, dentro de la condición.
  2. Refinar espacio

    Con h=1/20 y la misma Δt, λ=8/5.

    Para conservar 2/5 se reduce Δt a 1/1000, una cuarta parte.

  3. Implícito

    Con Euler implícito y λ=2, la alternancia tiene factor 1/9, estable. Un paso grande aún puede ser inadecuado para la precisión temporal solicitada.

Resultado

La elección espacial y la temporal están acopladas en el esquema explícito.

Lo comprobamos

La fórmula estudiada tiene orden temporal uno y espacial dos bajo regularidad suficiente.

8. El refinamiento útil controla el error de resolución y el redondeo

El error total reúne discretización, aproximación del sistema, evaluación de datos y aritmética. Si una fuente deja de ser pequeña respecto de otra, los cocientes entre mallas dejan de mostrar el orden formal. Se registran norma, residual original, tolerancias y tipo de frontera.

Un modelo ilustrativo de compromiso es E(h)=h²+ε/h²: el primer término baja al refinar, el segundo crece. No es una ley universal para cada algoritmo; explica por qué conviene medir en lugar de afirmar que el paso más pequeño siempre será mejor.

Una tolerancia lineal puede dominar a la discretización

Para el ejemplo de cuarto grado, N par, h=0,01; residual del sistema LₕU=f acotado por 0,001.

  1. Discretización

    El error nodal de la resolución exacta discreta es h²/4=0,000025.
  2. Resolución

    La cota del inverso permite hasta 0,001/8=0,000125 de error adicional por resolver aproximadamente.

  3. Ajustar

    Esa cota adicional es cinco veces la discretización. Con residual 10⁻⁶, baja a 0,000000125; entonces el refinamiento se estudia con menos interferencia de esa fuente.

Resultado

Una tolerancia fija puede ocultar el orden que buscamos medir.

Lo comprobamos

Estas son cotas de componentes del error; no se afirma que siempre alcancen simultáneamente sus máximos.

Un modelo de redondeo tiene un punto de equilibrio

Usar E(h)=h²+10⁻⁸/h², h>0.

  1. Paso 0,1

    E=0,01+0,000001=0,010001.
  2. Paso 0,01

    E=0,0001+0,0001=0,0002.
  3. Paso 0,001

    E=0,000001+0,01=0,010001.

    El mínimo de este modelo se alcanza en h=0,01, cuando ambos términos se igualan.

Resultado

Seguir refinando puede aumentar un error combinado.

Lo comprobamos

El equilibrio depende de las constantes y del modelo de error de cada cálculo.

Ahora lo resolvemos de principio a fin.

EJEMPLO RESUELTO 1

Un solo nodo interior

Resuelve −u″=2 en [0,1], u(0)=u(1)=0 con h=1/2.

  1. Malla

    Con h=1/2, los nodos son cero, un medio y uno. Las dos fronteras valen cero; solo U₁ es incógnita.

  2. Operador

    Para −u″=2, la fila es (2U₁−U₀−U₂)/h²=2.
  3. Datos conocidos

    Sustituir fronteras cero y h²=1/4 da 2U₁/(1/4)=2.
  4. Resolver

    8U₁=2, luego U₁=1/4.
  5. Control continuo

    u=x(1−x) cumple −u″=2 y ambos extremos cero; en x=1/2 vale 1/4.
  6. Alcance

    La fórmula centrada es exacta sobre esa cuadrática en los nodos. Unir solo esos puntos por rectas no vuelve exacta toda la parábola entre nodos.

Resultado

u₁=1/4.

Lo comprobamos

La solución exacta u=x(1−x) también vale 1/4 en el centro.

EJEMPLO RESUELTO 2

Aproximar una curvatura

Para u(x)=x², x=1, h=1/2, calcula la segunda diferencia.

  1. Fórmula

    La segunda diferencia es [u(x−h)−2u(x)+u(x+h)]/h².

  2. Argumentos

    Con x=1 y h=1/2 se evalúa en 1/2, uno y 3/2.
  3. Valores

    Para u=x² son 1/4, uno y 9/4.
  4. Numerador

    1/4−2+9/4=1/2.
  5. Dividir

    h²=1/4; (1/2)/(1/4)=2.
  6. Control

    La derivada exacta u″ es dos. La identidad (x−h)²−2x²+(x+h)²=2h² explica por qué esta fórmula es exacta para cualquier h no nula sobre una cuadrática.

Resultado

Aproximación 2.

Lo comprobamos

La derivada exacta u″ es 2.

Los símbolos que hemos usado

Malla / h
Conjunto de puntos de cálculo / distancia uniforme entre vecinos.
uᵢ / Uᵢ
Valor exacto u(xᵢ) / valor numérico calculado.
Stencil
Conjunto de vecinos y coeficientes que aproxima un operador.
Error de truncamiento
Defecto al introducir la solución exacta en el operador discreto.
Consistencia
El operador discreto aproxima el continuo al refinar.
Estabilidad
Control uniforme de cómo los defectos o perturbaciones afectan a la solución discreta.

Ahora te toca a ti.

Tienes cinco ejercicios. Inténtalos; las pistas te ayudan si te atascas.

Practicar este tema →

Para ampliar