Práctica 1. Modelación, método de Euler, error y convergencia

Curso: Modelación Matemática y Computacional I HPC

Tema: Solución numérica de problemas de valor inicial mediante el método de Euler

Pregunta central

¿Cómo se aproxima un problema de valor inicial mediante el método de Euler y cómo se estudian su error y convergencia?

Apuntes rápidos

Material de clase y actividades

1. Problemas de valor inicial

Muchas leyes de evolución pueden expresarse mediante una ecuación diferencial ordinaria de primer orden:

La ecuación describe la razón de cambio de una cantidad con respecto al tiempo . Sin embargo, esta ecuación por sí sola generalmente representa una familia de soluciones. Para seleccionar una solución particular necesitamos conocer el valor de la variable en algún instante:

El problema

se denomina problema de valor inicial (PVI).

En un método numérico no intentamos construir de inmediato una fórmula válida para todo . En su lugar, aproximamos la solución en una colección de puntos

Si utilizamos una malla uniforme,

donde se llama tamaño de paso.

La notación

representará una aproximación numérica de la solución exacta . En general,

Esta diferencia será central en el estudio del error.

2. Método de Euler

Supongamos, por un momento, que la solución exacta es suficientemente suave. Queremos relacionar el valor de la solución en con su valor un paso después,

La expansión de Taylor de alrededor de es

Como ,

Pero el PVI nos dice que

Por tanto,

y podemos escribir

Si conservamos solamente los dos primeros términos,

Finalmente, en el cálculo numérico sustituimos los valores exactos desconocidos por aproximaciones:

Esta es la fórmula del método de Euler explícito.

2.1 Interpretación geométrica

En el punto , la ecuación diferencial proporciona una pendiente

La recta con esa pendiente que pasa por es

Si evaluamos esta recta en ,

de donde

Actividad 1. Deducción

Sin copiar simplemente la fórmula final:

  1. Escribe la expansión de Taylor de alrededor de hasta el término de segundo orden.
  2. Sustituye utilizando la ecuación diferencial.
  3. Identifica el término que se omite al construir Euler.
  4. Explica con tus palabras por qué Euler puede interpretarse como una aproximación lineal local.
  5. Indica qué información requiere Euler para calcular .

La respuesta debe mostrar el desarrollo, no únicamente .

Actividad 1

Escribe aquí tu deducción. Utiliza ecuaciones en Markdown/LaTeX.

Hint: tu desarrollo debe contener un término proporcional a antes de llegar a la fórmula de Euler.

3. Modelo: enfriamiento en un ambiente cuya temperatura cambia

Consideremos un recipiente con líquido caliente. Inicialmente,

El recipiente se encuentra en una habitación cuya temperatura no permanece constante. Durante las primeras cinco horas se aproxima mediante

donde:

  • se mide en horas;
  • se mide en grados Celsius.

La rapidez con la que cambia la temperatura del líquido es proporcional a la diferencia entre su temperatura y la temperatura del ambiente. La constante de intercambio térmico es

La ley de enfriamiento de Newton puede expresarse como

¿Por qué aparece un signo negativo? Si el líquido está más caliente que el ambiente, entonces

y físicamente esperamos que la temperatura del líquido disminuya:

El signo negativo garantiza precisamente este comportamiento.

3.1 Construcción del PVI

Sustituyendo

obtenemos

con condición inicial

Por tanto, nuestra función es

Actividad 2. Interpretación del modelo

Responde antes de programar:

  1. ¿Cuál es la variable independiente?
  2. ¿Cuál es la variable dependiente?
  3. ¿Cuáles son las unidades de ?
  4. Verifica dimensionalmente que tiene las mismas unidades que .
  5. ¿Qué ocurriría con el signo de si el líquido estuviera más frío que el ambiente?
  6. ¿Qué representa físicamente el caso en un instante determinado?

Actividad 2

Hint: Si en un instante obliga a que ambas temperaturas permanezcan iguales posteriormente cuando cambia con .

4. Aplicación del método de Euler

Antes de escribir código, realicemos algunos pasos manualmente.

Para nuestro modelo,

Utilizaremos

Actividad 3. Iteraciones

Calcula manualmente usando .

Completa una tabla de la forma:

0090
1
2
3

Muestra al menos una sustitución completa en la fórmula de Euler.

Actividad 3

Completa la tabla anterior. Verifica tus cuatro iteraciones antes de compararlas con el código.

Incluye aquí al menos una sustitución completa:

5. Algoritmo

Un método numérico debe traducirse en una secuencia inequívoca de operaciones.

Para una función general

Euler requiere repetir:

Un esquema básico sería:

  1. guardar y ;
  2. evaluar ;
  3. calcular ;
  4. avanzar el tiempo;
  5. guardar los nuevos valores;
  6. repetir hasta llegar a .

Si

no es entero, no podemos efectuar exactamente iteraciones. En esta práctica escogeremos valores de que dividan exactamente el intervalo, pero en una implementación general este detalle debe tratarse explícitamente.

Actividad 4. Pseudocódigo

Escribe pseudocódigo para una función que reciba

y regrese dos arreglos con los tiempos y las aproximaciones numéricas.

Tu algoritmo debe indicar dónde:

  • se calcula el número de pasos;
  • se crea la malla;
  • se almacena la condición inicial;
  • se evalúa la pendiente;
  • se actualiza ;
  • se regresan los resultados.

Actividad 4

función Euler(f, t0, y0, tf, h)
 
    # Completar
 
fin función

No escribas todavía Python: primero representa la lógica del algoritmo.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

6. Implementación del modelo

La función matemática del problema es

Separar el modelo del método numérico es una buena práctica de programación:

  • f(t,T) describe el fenómeno;
  • euler(...) describe el método numérico.

Así, la misma implementación de Euler puede utilizarse posteriormente con otras ecuaciones diferenciales.

def f_temperatura(t, T):
    # TODO 1:
    # Calcular la temperatura ambiente Ta(t) = 20 + 2t.
    Ta = None
 
    # TODO 2:
    # Regresar -0.6*(T - Ta).
    return None
 
 
def euler(f, t0, y0, tf, h):
    """Implementa aquí el método de Euler explícito."""
 
    # TODO 1: calcular cuántos pasos se necesitan.
    N = None
 
    # TODO 2: crear arreglos para t y y con N+1 elementos.
    # t = ...
    # y = ...
 
    # TODO 3: almacenar la condición inicial.
    # t[0] = ...
    # y[0] = ...
 
    # TODO 4: construir un ciclo for.
    # En cada paso:
    #   a) calcular t[n+1]
    #   b) evaluar f(t[n], y[n])
    #   c) calcular y[n+1]
 
    # TODO 5: regresar ambos arreglos.
    # return t, y
 
    raise NotImplementedError("Completa la función euler antes de continuar.")

Comprobació

Cuando termines tu función, resuelve solamente hasta con y compara los primeros cuatro pasos con los resultados que calculaste manualmente.

Si no coinciden, corrige el código antes de continuar.

# Descomenta cuando tu función euler esté terminada.
 
# t_test, T_test = euler(f_temperatura, 0.0, 90.0, 2.0, 0.5)
# pd.DataFrame({"t": t_test, "T_Euler": T_test})

7. Primer experimento: efecto del tamaño de paso

Resolveremos el modelo en

con diferentes tamaños de paso:

Antes de ejecutar el experimento, piensa:

  • al disminuir , ¿aumentará o disminuirá el número de pasos?
  • ¿esperas que las curvas numéricas se acerquen entre sí?
  • ¿por qué una malla más fina suele producir una mejor aproximación?
h_values = [1.0, 0.5, 0.25, 0.125]
 
# TODO:
# 1. crea una figura;
# 2. recorre h_values;
# 3. llama a tu función euler;
# 4. grafica cada aproximación;
# 5. agrega título, etiquetas, leyenda y cuadrícula.
 
# Tu código aquí.

Análisis

Explica:

  1. qué ocurre con el número de pasos cuando disminuye;
  2. qué cambia en las curvas;
  3. por qué este experimento sugiere pero todavía no demuestra cuantitativamente convergencia;
  4. cuál es el costo computacional de disminuir .

8. Solución analítica del modelo

Para estudiar el error necesitamos una referencia. En este caso particular el PVI puede resolverse analíticamente.

Partimos de

Distribuyendo,

por lo que

Esta es una ecuación diferencial lineal de primer orden

con

El factor integrante es

Multiplicando toda la ecuación por ,

El lado izquierdo es

Por tanto,

Integramos ambos lados:

Al efectuar la integración y aplicar , se obtiene

Verificación de la solución exacta

Una fórmula obtenida analíticamente no debe aceptarse sin revisión.

Condición inicial

Derivada

Al sustituir en

se obtiene la misma expresión. Por tanto, la función satisface tanto la ecuación diferencial como la condición inicial.

def solucion_exacta(t):
    """Solución analítica del modelo de temperatura."""
    return 56/3 + 2*t + (214/3)*np.exp(-0.6*t)
 
print("T_exacta(0) =", solucion_exacta(0.0))
print("T_exacta(5) =", solucion_exacta(5.0))

9. Error

Es importante distinguir conceptos.

9.1 Error de la aproximación en un punto

Si es el valor exacto y es la aproximación numérica,

El error absoluto es

El error relativo es

siempre que .

El error porcentual es

Estos errores comparan directamente la solución numérica con una referencia exacta. No deben confundirse automáticamente con el error local de truncamiento.

10. Error local de truncamiento

Con Taylor:

para algún .

Euler conserva

y omite el término restante. El error cometido en un solo paso ideal, suponiendo que iniciamos ese paso desde el valor exacto, es entonces proporcional a

Por eso decimos que el error local de truncamiento de Euler es

El símbolo no significa que el error sea exactamente . Significa que, para suficientemente pequeño, su magnitud puede acotarse por una constante multiplicada por .

11. Error global y orden del método

En un intervalo de longitud fija,

el número de pasos es aproximadamente

Por tanto,

Aunque cada paso introduce un error local de orden , realizamos del orden de pasos. De manera intuitiva,

Bajo las hipótesis matemáticas apropiadas, el error global del método de Euler satisface

Por esta razón decimos que Euler es un método de primer orden.

# Construye una tabla de errores en t=5.
 
h_values_error = np.array([1.0, 0.5, 0.25, 0.125, 0.0625, 0.03125])
T5_exacta = solucion_exacta(5.0)
 
filas = []
 
for h in h_values_error:
    # TODO:
    # t_num, T_num = ...
    # error_abs = ...
    # error_rel = ...
    # error_porcentual = ...
    #
    # filas.append([...])
    pass
 
# TODO: crea un DataFrame con las columnas:
# h, N, T_Euler(5), T_exacta(5),
# Error absoluto, Error relativo, Error (%)
 
# tabla_error = ...
# tabla_error

Actividad 5. Diferencia entre error local y error global

Explica con tus propias palabras por qué no existe contradicción entre

para el error local y

para el error global.

Tu explicación debe distinguir explícitamente:

  • un solo paso;
  • muchos pasos;
  • acumulación/propagación del error;
  • intervalo de integración fijo.

Actividad 5

Desarrolla la explicación en un párrafo matemáticamente preciso. Deben aparecer las expresiones

No basta con multiplicarlas: explica qué representa cada una.

12. Convergencia

Decimos que un método numérico converge si, al refinar la malla,

la solución numérica se aproxima a la solución exacta:

La convergencia no significa simplemente que “la gráfica se ve bien”. Es una propiedad cuantitativa.

Si el error se comporta como

el exponente es el orden de convergencia.

Para dos pasos consecutivos y ,

Dividiendo,

Así,

Para Euler esperamos

13. Gráfica log-log y orden experimental

A partir de

tomamos logaritmos:

Esto tiene la forma de una recta

con

Por tanto, si graficamos el error contra en escala log-log, la pendiente debería aproximarse al orden del método.

Para Euler esperamos una pendiente cercana a

# ORDEN EXPERIMENTAL
#
# Usa los errores de tu tabla para calcular:
#
# p = ln(E(h)/E(h/2)) / ln(2)
#
# Construye una tabla con:
# h, E(h), E(h)/E(h/2), p experimental
 
# Tu código aquí.
# GRÁFICA LOG-LOG
#
# 1. Obtén log(h) y log(E).
# 2. Usa np.polyfit(..., ..., 1) para ajustar una recta.
# 3. La primera componente del resultado será la pendiente.
# 4. Grafica los datos y reporta el valor de p.
 
# Tu código aquí.

Análisis de convergencia

Responde:

  1. ¿Los valores de se acercan a ?
  2. ¿La pendiente log-log es compatible con el orden teórico de Euler?
  3. ¿Qué relación empírica entre y sugieren tus resultados?
  4. ¿Por qué los cocientes pueden no ser exactamente iguales a ?

14. Preguntas

Responde con argumentos, no únicamente con una palabra o fórmula.

  1. ¿Qué información utiliza Euler para avanzar de a ?
  2. ¿Por qué se dice que es un método explícito?
  3. ¿De dónde surge matemáticamente la fórmula de Euler?
  4. ¿Qué relación existe entre la curvatura de una solución y el error local?
  5. ¿Cuál es la diferencia entre y ?
  6. ¿Cuál es la diferencia entre error absoluto, error relativo y error local de truncamiento?
  7. ¿Por qué Euler tiene error local pero error global ?
  8. ¿Qué significa afirmar que un método converge?
  9. ¿Qué significa afirmar que Euler es de orden uno?
  10. Si reducimos a la mitad, ¿qué esperamos que ocurra aproximadamente con el error global?
  11. ¿Qué costo computacional tiene reducir ?
  12. ¿Por qué el análisis de estabilidad impide elegir únicamente por comodidad?
  13. En el modelo térmico, ¿qué significado físico tiene el signo de ?
  14. ¿Qué ventajas tiene separar en el código la función del modelo y la implementación del método?

Respuestas

Dudas

  • [ ]

Conceptos para extraer

Referencias mencionadas

  • Notebook de clase: Práctica1_MMyCHPC_Euler_clase.ipynb

Resumen después de clase

En 3–5 oraciones, ¿qué aprendí y cómo se conecta con clases anteriores?

Acciones adicionales

  • Completar las actividades y celdas de código durante la clase.
  • Resolver dudas pendientes.