Introducción a las Ecuaciones diferenciales

Introducción

Las ecuaciones diferenciales son relaciones matemáticas entre una función y sus derivadas. Por ejemplo:

\[ y' = ky \qquad \qquad k \in \mathbb{R} \]

El estudio de estas ecuaciones se motiva por su capacidad para modelizar fenómenos naturales que se expresan naturalmente como ecuaciones diferenciales

Ejemplo 1.1 Un ejemplo conocido es el disparo parabólico, que surge de aplicar las leyes de Newton al movimiento de un proyectil. Un cañón dispara una bala desde la altura \(h_0\) con una velocidad de \(30 m/s\) y un ángulo \(\alpha\). Se puede expresar la posición de bala de cañón a través del siguiente sistema, donde el vector \((x(t), y(t))\) describe la posición de la bala en cada instante de tiempo \(t\):

\[ \begin{cases} y'' = -9.81 \\ x'' = 0 \\ y'(0) = 30 \sin \alpha \\ x'(0) = 30 \cos \alpha\\ y(0) = h_0\\ x(0) = 0 \\ \end{cases} \]

Código

v0 = 30;
g = 9.81;

// Altura máxima absoluta posible (disparo vertical desde altura máxima)
h_max_absolute = 40 + (v0 * v0) / (2 * g);

// Alcance máximo absoluto (ángulo óptimo desde h0=50)
alcance_max_absolute = (() => {
  const h0_max = 50;
  const ang_opt = 45 * Math.PI / 180; // aproximación
  const vx = v0 * Math.cos(ang_opt);
  const vy = v0 * Math.sin(ang_opt);
  const disc = vy * vy + 2 * g * h0_max;
  const t = (vy + Math.sqrt(disc)) / g;
  return vx * t;
})();

// Cálculos
ang_rad = ang * Math.PI / 180;
vx = v0 * Math.cos(ang_rad);
vy = v0 * Math.sin(ang_rad);

t_hit = (() => {
  const a = -0.5 * g;
  const b = vy;
  const c = h0;
  const disc = b * b - 4 * a * c;
  if (disc <= 0) return 0;
  const r1 = (-b + Math.sqrt(disc)) / (2 * a);
  const r2 = (-b - Math.sqrt(disc)) / (2 * a);
  return Math.max(r1, r2);
})();

alcance = vx * t_hit;
h_max = h0 + (vy * vy) / (2 * g);
periodo = t_hit + 0.8;

// Trayectoria completa
trajectory = {
  const pts = [];
  const dt = Math.max(t_hit / 220, 0.01);
  for (let t = 0; t <= t_hit; t += dt) {
    const x = vx * t;
    const y = h0 + vy * t - 0.5 * g * t * t;
    pts.push({t, x, y: Math.max(0, y)});
  }
  return pts;
};

// Animación
t_anim = (now / 500) % periodo;
t_live = Math.min(t_anim, t_hit);
x_live = vx * t_live;
y_live = Math.max(0, h0 + vy * t_live - 0.5 * g * t_live * t_live);

Plot.plot({
  width: 700,
  height: 300,
  marginLeft: 60,
  marginBottom: 55,
  grid: true,
  x: {
    label: "Distancia horizontal (m) →",
    domain: [-5, alcance_max_absolute * 1.15]
  },
  y: {
    label: "↑ Altura (m)",
    domain: [-2, h_max_absolute * 1.1]
  },
  marks: [
    // Suelo
    Plot.ruleY([0], {stroke: "#654321", strokeWidth: 4}),
    
    // Trayectoria parabólica
    Plot.line(trajectory, {
      x: "x",
      y: "y",
      stroke: "#17a2b8",
      strokeWidth: 2.5,
      strokeDasharray: "5,3"
    }),
    
    // Proyectil en movimiento
    Plot.dot([{x: x_live, y: y_live}], {
      x: "x",
      y: "y",
      fill: "#e74c3c",
      stroke: "#c0392b",
      strokeWidth: 2,
      r: 7
    }),
    
    // Línea de alcance
    Plot.ruleX([alcance], {
      stroke: "#27ae60",
      strokeWidth: 2,
      strokeDasharray: "6,4"
    }),
    
    // Texto con información
    Plot.text([{x: alcance / 2, y: h_max * 0.95}], {
      text: [`Alcance: ${alcance.toFixed(1)} m | Altura máx: ${h_max.toFixed(1)} m`],
      fill: "#2c3e50",
      fontSize: 13,
      fontWeight: "bold"
    })
  ]
})

Observaciones:

Alcance horizontal: \(R \approx\) m
Tiempo de vuelo: \(T \approx\) s
Altura máxima: \(h_{\max} \approx\) m

Este problema conocido por el lector introduce algunas preguntas clave de esta asignatura:

¿El problema tiene solución? ¿Es única?
¿Se puede resolver de forma explícita?
¿Se puede aproximar numéricamente?
¿Como varían las soluciones en función de los parámetros como gravedad y altura inicial?

Nuestra capacidad de entender y predecir fenómenos del mundo real depende directamente de comprender y analizar las ecuaciones diferenciales que lo describen.

Ejercicio 1.1 Resolver separadamente las ecuaciones del problema anterior \(y'' = -9.81\), \(x'' = 0\) y entender como se relacionan con el problema

Los tres enfoques distintos a las ecuaciones diferenciales

Las preguntas del ejemplo anterior adelantan los tres enfoques de estudio principales:

Analítico: centrado en calcular una expresión exacta de las soluciones
Cualitativo: centrado en describir el comportamiento de las soluciones sin resolverlas
Numérico: centrado en aproximar soluciones controlando el error

Estos apuntes se centran en métodos de resolución analítica y propiedades de ecuaciones particulares. Cursos posteriores profundizarán en análisis cualitativo (teoría de sistemas dinámicos), métodos numéricos y ecuaciones en derivadas parciales.

Antes de desarrollar la teoría formal, veamos qué tipo de preguntas se pueden responder sobre una ecuación diferencial simple.

Ejemplo informal: crecimiento bacteriano

A continuación se realiza un breve estudio de la ecuación \(y' = a y\). El ejemplo está diseñado para que el alumno intente resolver todos los apartados.

Ejemplo 1.2 La ecuación \(y' = a y\) es una ecuación lineal de primer orden homogenea. Será objeto de estudio a lo largo de la asignatura.

El objetivo de este ejemplo es presentar de manera intuitiva como se pueden estudiar estos problemas.

Sin tener la solución, ¿Cuando la solución es creciente/decreciente?

Solución 1.2. Aunque no se tiene la solución, la ecuación diferencial aporta información sobre su derivada. En particular la función es creciente si y solo si \(y'(x) > 0\).

\[ \boxed{\;y'(x) > 0 \quad \Longleftrightarrow \quad a y(x) > 0 \quad \Longleftrightarrow \quad \operatorname{sgn}(y)= \operatorname{sgn}(a)\;} \]

Esto quiere decir que depende de los valores de a: - Si \(a > 0\), las soluciones positivas son crecientes y las negativas decrecientes. - Si \(a < 0\), las soluciones positivas son decrecientes y las negativas crecientes.

A pesar de que no se conoce la solución, se pueden deducir algunas de sus caracteristicas a partir de la ecuación:

Comprobar que \(y(x) = Ce^{a x}\) con \(C \in \mathbb{R}\) es solución de \(y' = a y\)

Solución 1.2. La función \(y(x) = Ce^{a x}\) es solución de \(y' = a y\). Para demostrarlo se sustituye en la expresión: \[ \boxed{\frac{d}{dx} (Ce^{a x}) = a Ce^{a x} = a y} \]

Sin embargo, surge una pregunta: ¿Existen otras soluciones?

Demostrar que todas las soluciones son de la forma \(y(x) = Ce^{a x}\) con \(C \in \mathbb{R}\). (Pista: multiplicar por \(e^{-ax}\))

Solución 1.3. Si multiplicamos la ecuación por la función \(e^{-a x}\): \[ \boxed{e^{-a x} y' - a e^{-a x} y = 0 \implies \frac{d}{dx}(e^{-a x} y) = 0} \]

Por lo tanto, \(e^{-a x} y\) es constante y todas las soluciones son de la forma \(y(x) = Ce^{a x}\).

La función \(e^{-a x}\) es el factor integrante de la ecuación que proporciona un método de resolución que se estudiará más adelante

La siguiente pregunta es: ¿Qué caracteriza esa constante?

Encontrar el valor de la función que caracteriza la constante.

Solución 1.4. Como en el ejemplo del disparo, se puede tener un punto de la solución como condición de inicio: \(y(0) = y_0\). Entonces la solución es: \[ \boxed{y_0 = y(0) = C e^{a 0} = C \quad \Longrightarrow \quad y(x) = y_0 e^{a x}} \]

Los resultados obtenidos hasta ahora nos ayudan a intuir que las ecuaciones tienen como solución familias de funciones y que las condiciones que se imponen determinan subconjuntos particulares.

La siguiente pregunta es: ¿Cómo varían las soluciones respecto de a?

Describir como varían las soluciones respecto al parámetro a.

Solución 1.5. El comportamiento depende del parámetro \(a\):

Si \(a > 0\), las soluciones \(y(x) = y_0 e^{a x}\) son exponenciales crecientes.
Si \(a < 0\), las soluciones \(y(x) = y_0 e^{a x}\) son exponenciales decrecientes.
Si \(a = 0\), las soluciones \(y(x) = y_0\) son constantes.

Nota: el crecimiento se podía deducir directamente en el primer apartado.

Para estudiar el espacio de soluciones es útil representar todas las soluciones en el plano \((x, y)\). Este espacio se denomina plano de fases.

Representar todas las soluciones en el plano de fases.

Solución 1.6. En este caso, las soluciones son curvas exponenciales:

Código

viewof a_fase = Inputs.range([-2, 2], {
  value: 0.5,
  step: 0.01,
  label: "a"
})

// Generar muchas condiciones iniciales
condiciones_iniciales_fase = Array.from({length: 21}, (_, i) => -10 + i);

// Generar órbitas para cada condición inicial
orbitas_fase = condiciones_iniciales_fase.map(y0 => {
  const puntos = [];
  const tmax = 8;
  const dt = 0.05;
  for (let t = 0; t <= tmax; t += dt) {
    const y = y0 * Math.exp(a_fase * t);
    puntos.push({t, y});
  }
  return {y0, puntos};
});

// Gráfico de plano de fases (t, y)
Plot.plot({
  width: 700,
  height: 400,
  marginLeft: 60,
  marginBottom: 55,
  grid: true,
  x: {
    label: "t →",
    domain: [0, 8]
  },
  y: {
    label: "↑ y(t)",
    domain: [-45, 45]
  },
  marks: [
    Plot.ruleY([0], {stroke: "#aaa", strokeWidth: 1.5}),
    Plot.ruleX([0], {stroke: "#aaa", strokeWidth: 1.5}),
    ...orbitas_fase.map(orb => Plot.line(orb.puntos, {
      x: "t",
      y: "y",
      stroke: orb.y0 === 0 ? "#27ae60" : (a_fase > 0 ? "#c2300fff" : a_fase < 0 ? "#27ae60" : "#888"),
      strokeWidth: orb.y0 === 0 ? 3.5 : 2.2,
      opacity: orb.y0 === 0 ? 1 : 0.7
    })),
    Plot.dot(orbitas_fase.map(orb => ({t: 0, y: orb.y0})), {
      x: "t",
      y: "y",
      fill: "#222",
      r: 4,
      opacity: 0.7
    })
  ],
  style: {
    background: "#fafafa"
  }
})

En el plano se observa que hay una solución constante \(y=0\). Estas soluciones se llaman equilibrios o soluciones estacionarias.

Encontrar todos los equilibrios o soluciones constantes.

Solución 1.7. Soluciones particulares y de interés son las soluciones constates, llamadas puntos de equilibrio. Una solución es constante \(y(x) = k \in \mathbb{R}\) si y solo si \(y'(x) = 0\). Por lo tanto las soluciones de: \[ y' = 0 = a y \quad \Longleftrightarrow \quad a = 0 \lor y = 0 \]

Si \(a \neq 0\), el único equilibrio es la función \(0\).
Si \(a = 0\), todas las soluciones son equilibrios.

La interpretación geómetrica de la derivada es la pendiente de la recta tangente a la solución. Esto quiere decir que en cada punto la ecuación define la tangente a una o varias soluciones.

Esto nos permite dibujar un campo vectorial asociado en el que en cada punto el vector tiene la pendiente definida por la ecuación.

Definir el campo vectorial asociado a la ecuación diferencial y dibujarlo.

Solución 1.8. El campo de direcciones muestra la pendiente \(y' = a y\) en cada punto \((x, y)\). Este permitirá visualizar el comportamiento local de las soluciones.

Código

viewof a = Inputs.range([-3, 3], {
  value: -1,
  step: 0.1,
  label: "a"
})
{
  // Dominio
  let x_min = -3, x_max = 3;
  let y_min = -3, y_max = 3;

  // Campo
  let densidad = 20;
  let longitud = 0.35;
  let normalizar = true;

  // Sistema: x' = 1, y' = a*y
  let f = (x, y, t) => 1;
  let g = (x, y, t) => a * y;

  // Color por estabilidad (tiempo hacia delante)
  let colorCampo = (a < 0) ? "green" : (a > 0) ? "red" : "#777";

  // Campo vectorial (2D)
  let calcular_campo_vectorial = (f, g, x_min, x_max, y_min, y_max, densidad, longitud, normalizar) => {
    let vectores = [];
    let dx = (x_max - x_min) / densidad;
    let dy = (y_max - y_min) / densidad;

    for (let x = x_min; x <= x_max; x += dx) {
      for (let y = y_min; y <= y_max; y += dy) {
        let vx = f(x, y, 0);
        let vy = g(x, y, 0);
        if (!isFinite(vx) || !isFinite(vy)) continue;

        let mag = Math.hypot(vx, vy);
        if (mag === 0) continue;

        let scale = normalizar
          ? Math.min(dx, dy) * longitud
          : Math.min(dx, dy) * longitud * Math.min(1, mag);

        let dx_vec = (vx / mag) * scale;
        let dy_vec = (vy / mag) * scale;

        vectores.push({ x1: x, y1: y, x2: x + dx_vec, y2: y + dy_vec });
      }
    }
    return vectores;
  };

  let campo = calcular_campo_vectorial(
    f, g,
    x_min, x_max,
    y_min, y_max,
    densidad,
    longitud,
    normalizar
  );

  return Plot.plot({
    width: 700,
    height: 400,
    marginLeft: 60,
    marginBottom: 60,
    grid: true,
    x: { label: "x →", domain: [x_min, x_max] },
    y: { label: "↑ y", domain: [y_min, y_max] },
    marks: [
      Plot.ruleX([0], { stroke: "black", strokeWidth: 1.5, opacity: 0.25 }),
      Plot.ruleY([0], { stroke: "black", strokeWidth: 1.5, opacity: 0.25 }),

      Plot.arrow(campo, {
        x1: "x1", y1: "y1",
        x2: "x2", y2: "y2",
        stroke: colorCampo,
        strokeWidth: 1.8,
        headLength: 8,
        headAngle: 25,
        opacity: (a === 0) ? 0.35 : 0.75
      })
    ]
  });
}

Para definir el campo es necesario un vector \((\alpha, \beta)\) que tenga como pendiente \(ay\).

Como la pendiente es \(\dfrac{\beta}{\alpha}\), el vector más simple es \((1, ay)\).

Interpretación: Cada flecha indica la dirección y pendiente de la solución en ese punto.

El dibujo se ha hecho numéricamente evaluando en todos los puntos aunque se verán tecnicas para poder aproximarlo sin la necesidad de evaluar tantos puntos.

Además se aprecia que si \(a >0\) las soluciones se alejan del \(0\) mientras que si es negativo se acercan. Esta propriedad de la solución \(0\) se estudiará como estabilidad

Formalizar la intuición anterior usando la expresión de la solución. ¿Cuando convergen las soluciones a 0?

Solución 1.10. Como las soluciones son exponenciales. El límite a infinito converge a \(0\) si el argumento es negativo y convergen a \(\pm \infty\) si es positivo.

Formalmente se dirá que:

Para \(a < 0\), el equilibrio \(y = 0\) es estable (las soluciones tienden a cero).
Para \(a > 0\), el equilibrio \(y = 0\) es inestable (las soluciones se alejan de cero).

En teoría de sistemas dinámicos estos cambios de comportamiento son muy importantes. Un cambio significativo en el comportamiento de las soluciones puede hacer que un modelo deje de ser representativo o un aviso de que el fenómeno modelado cambiará su comportamiento. Estos eventos se denominan bifurcaciones.

Además el campo de direcciones permite aproximar numéricamente las soluciones. Si aproximamos la solución en un punto por la tangente a ese punto se obtiene el método de Euler.

Para ello, dividimos el espacio en un conjunto de \(N\) puntos equidistantes y de forma iterada aproximamos en cada punto la función por la recta tangente definida por la EDO.

Formalizar la intuición anterior y aproximar \(y' = y\) en \([0,1]\) donde \(y(0)=1\).

Solución 1.10. Para aproximar la solución del problema \(y' = y\), \(y(0) = 1\) en el intervalo \([0,1]\) mediante el método de Euler, dividimos el intervalo en \(N\) pasos de tamaño \(h = 1/N\).

Además para aproximar un punto se usa el punto anterior y la recta tangente a definida por la EDO \(y_{n+1} = y_{n} + h \dot y_n = y_{n-1} + h y_{n-1}\) \[ \begin{cases} \displaystyle x_n = n h = \frac{n}{N} \\ \displaystyle y_0 = 1 \\ \displaystyle y_n = y_{n-1} + h y_{n-1} = (1 + h) y_{n-1} \end{cases} \]

para \(n = 1, 2, \ldots, N\).

Por inducción se obtiene la expresión explícita: \[ y_n = (1 + h)^n = \left(1 + \frac{1}{N}\right)^n \]

Código

viewof N = Inputs.range([1, 25], {
  value: 8,
  step: 1,
  label: "N"
})

viewof dibujarPoligonal = Inputs.toggle({
  value: true,
  label: "Dibujar poligonal"
})

viewof mostrarExacta = Inputs.toggle({
  value: false,
  label: "Mostrar solución analítica (e^x)"
})

viewof animar = Inputs.toggle({
  value: true,
  label: "Animar"
})

// “reloj” como generador asíncrono (compatible con Quarto OJS)
frame = animar
  ? (async function* () {
      let m = 0;
      const total = 4 * N;   // 4 fases por paso
      while (true) {
        yield m;
        m = (m + 1) % (total + 1);
        await new Promise((r) => setTimeout(r, 250));
      }
    })()
  : 4 * N;

{
  const x0 = 0, x1 = 1;
  const y0 = 1;
  const h = (x1 - x0) / N;

  // Euler explícito: y_{k+1} = y_k + h*y_k
  const datos = [];
  let x = x0, y = y0;
  datos.push({ x, y, k: 0 });

  for (let k = 0; k < N; k++) {
    y = y + h * y;
    x = x + h;
    datos.push({ x, y, k: k + 1 });
  }

  // Exacta (opcional)
  const exacta = Array.from({ length: 401 }, (_, i) => {
    const xx = x0 + (x1 - x0) * (i / 400);
    return { x: xx, y: Math.exp(xx) };
  });

  // --- Fases de animación ---
  // frame m -> k y fase:
  // fase 0: solo punto actual
  // fase 1: aparece segmento (recta de Euler)
  // fase 2: aparece nuevo punto (segmento aún visible)
  // fase 3: se borra segmento, quedan puntos (listo para repetir)
  const m = animar ? frame : 4 * N;

  const k = Math.min(N, Math.floor(m / 4));
  const fase = m % 4;

  // cuántos puntos se ven
  let ultimoIndiceVisible;
  if (!animar) {
    ultimoIndiceVisible = N;
  } else {
    if (k >= N) {
      ultimoIndiceVisible = N;
    } else {
      // en fase 0 y 1: hasta k; en fase 2 y 3: hasta k+1
      ultimoIndiceVisible = (fase <= 1) ? k : (k + 1);
    }
  }

  const visibles = datos.slice(0, ultimoIndiceVisible + 1);

  // segmento del paso k (solo en fase 1 y 2)
  let segmento = [];
  if (animar && k < N && (fase === 1 || fase === 2)) {
    const pk = datos[k];
    const x1s = pk.x + h;
    const y1s = pk.y + h * pk.y;
    segmento = [{ x: pk.x, y: pk.y }, { x: x1s, y: y1s }];
  }

  // punto actual resaltado (fase 0/1)
  let puntoActual = [];
  if (animar && k <= N && (fase === 0 || fase === 1)) {
    puntoActual = [datos[k]];
  }

  // punto nuevo resaltado (fase 2/3)
  let puntoNuevo = [];
  if (animar && k + 1 <= N && (fase === 2 || fase === 3)) {
    puntoNuevo = [datos[k + 1]];
  }

  return Plot.plot({
    width: 700,
    height: 420,
    marginLeft: 60,
    marginBottom: 55,
    grid: true,
    x: { label: "x →", domain: [x0, x1] },
    y: { label: "↑ y", domain: [0.8, 2.9] },
    marks: [
      ...(mostrarExacta
        ? [Plot.line(exacta, { x: "x", y: "y", stroke: "#9ff383ff", strokeWidth: 2, opacity: 0.7 })]
        : []),

      // poligonal (solo con lo ya “aceptado/visible”)
      ...(dibujarPoligonal
        ? [Plot.line(visibles, { x: "x", y: "y", stroke: "black", strokeWidth: 2, opacity: 0.9 })]
        : []),

      // segmento que aparece y luego se borra
      ...(segmento.length
        ? [Plot.line(segmento, { x: "x", y: "y", stroke: "steelblue", strokeWidth: 3, opacity: 0.95 })]
        : []),

      // puntos visibles
      Plot.dot(visibles, { x: "x", y: "y", r: 3, fill: "black" }),

      // resaltados
      ...(puntoActual.length
        ? [Plot.dot(puntoActual, { x: "x", y: "y", r: 6, fill: "steelblue" })]
        : []),
      ...(puntoNuevo.length
        ? [Plot.dot(puntoNuevo, { x: "x", y: "y", r: 6, fill: "steelblue" })]
        : [])
    ]
  });
}

En particular, para \(x = 1\) (\(n = N\)): \[ y_N = \left(1 + \frac{1}{N}\right)^N \]

Comparando con la solución exacta \(y(1) = e^1 = e\), se observa que \(y_N\) converge a \(e\).

La ecuación como modelo

En este caso la ecuación modela dos fénomenos distintos para valores de \(a\) positivos o negativos.

Para valores de \(a\) positivos, es modelo muy sencillo para describir el crecimiento de una población de bacterias.

El crecimiento es proporcional a la cantidad de bacterias \(y' = ay\) (si se duplica la población la velocidad también duplica).

Pregunta: ¿Qué limitaciones podría tener este modelo?

Para valores negativos, modeliza decricimiento. Por ejemplo se usa para medir la cantidad de un elemento radioactivo.

El decadimento nuclear se comporta exactamente como esta ecuación

Notación y clasificación de ecuaciones

Notación de derivadas

Existen varias notaciones comunes para expresar derivadas en ecuaciones diferenciales:

Notación de Lagrange: \[y', \quad y'', \quad y'''\]

Notación de Newton: \[\dot{y}, \quad \ddot{y}, \quad \dddot{y}\]

Notación de Leibniz: \[\frac{dy}{dt}, \quad \frac{d^2y}{dt^2}, \quad \frac{d^3y}{dt^3}\]

A lo largo del texto se utilizarán principalmente la notación de Lagrange y de Leibniz.

Además a lo largo del libro se usarán distintas letras para la variable independiente (t, x, s, etc.) y la función desconocida (y, u, v, etc.) dependiendo del contexto.

Clasificación de ecuaciones diferenciales

Las ecuaciones diferenciales se clasifican según el número de variables independientes:

Ecuaciones Diferenciales Ordinarias (EDO) La función desconocida depende de una única variable independiente. Las derivadas son ordinarias.

\[y' = ky, \quad y'' + \omega^2 y = 0\]

Ecuaciones en Derivadas Parciales (EDP) La función desconocida depende de dos o más variables independientes. Aparecen derivadas parciales.

\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}, \quad \nabla^2 u = 0\]

Definición 1.1 (Orden) El orden de una ecuación diferencial es el de la derivada de mayor orden que aparece en ella.

\[ y^{(n)} = f(t, y, y', y'', \ldots, y^{(n-1)}) \]

La ecuación anterior es de orden \(n\).

Ejemplo 1.3

\(y' = 3y^3\) es de orden 1
\(y'' + 2y' + y = 0\) es de orden 2
\(y''' - y' = \sin(t)\) es de orden 3
\(\frac{d^4y}{dt^4} + 2\frac{d^2y}{dt^2} = t^2\) es de orden 4

Definición 1.2 (Ecuación lineal) Una ecuación diferencial ordinaria de orden \(n\) es lineal si puede escribirse en la forma:

\[ a_n(t)y^{(n)} + a_{n-1}(t)y^{(n-1)} + \cdots + a_1(t)y' + a_0(t)y = g(t) \]

donde \(a_i(t)\) y \(g(t)\) son funciones que dependen solo de \(t\) (no de \(y\) ni sus derivadas).

Características clave:

La variable dependiente \(y\) y sus derivadas aparecen de forma lineal (no hay productos, potencias, ni composiciones)
Si \(g(t) = 0\), la ecuación es homogénea
Si \(g(t) \neq 0\), la ecuación es no homogénea

Ejemplo 1.4 Ecuaciones lineales:

\(y' + 2y = e^t\) (lineal de orden 1, no homogénea)
\(y'' - 3y' + 2y = 0\) (lineal de orden 2, homogénea)
\(t^2y'' + ty' - y = \ln t\) (lineal de orden 2, no homogénea)

Ecuaciones no lineales:

\(y' = y^2\)
\(y'' + (y')^2 + y = 0\)
\(y'' + \sin(y) = 0\)

Una ecuación diferencial ordinaria es autónoma si la variable independiente no aparece explícitamente en la ecuación.

Forma general de una EDO autónoma de orden 1: \[y' = f(y)\]

Forma general de una EDO autónoma de orden 2: \[y'' = f(y, y')\]

Ejemplo 1.5 Ecuaciones autónomas:

\(y' = ky\)
\(y' = y(1 - y)\)
\(y'' + \sin(y) = 0\)

Ecuaciones no autónomas:

\(y' = t + y\)
\(y'' + ty' + y = 0\)
\(y' = e^t y\)

Forma estándar y forma normal

Una ecuación diferencial ordinaria de orden \(n\) se puede expresar en forma estándar o forma normal:

Definición 1.3 (Forma estándar) \[ F\left(t, y, y', y'', \ldots, y^{(n)}\right) = 0 \]

La ecuación admite forma normal (o forma explícita) si es posible despejar la derivada de mayor orden:

Definición 1.4 (Forma Normal) \[ y^{(n)} = f\left(t, y, y', y'', \ldots, y^{(n-1)}\right) \]

Ejemplo 1.6 La ecuación \(y'' + 2y' + y = e^t\) en forma normal es: \[ y'' = e^t - 2y' - y \]

En forma estándar es:

\[ y'' - e^t + 2y' + y = 0 \]

Soluciones de ecuaciones

Definición 1.5 (Solución) Una función \(y: I \subset \mathbb{R} \to \mathbb{R}\) (o \(y: I \to \mathbb{R}^m\)) (derivable tantas veces como lo exija la ecuación) se dice solución de una EDO en el intervalo \(I\) si, al sustituir \(y(t)\) y sus derivadas en la ecuación, esta se verifica para todo \(t \in I\).

En particular, si la EDO está en forma normal \[ y^{(n)} = f\left(t, y, y', \ldots, y^{(n-1)}\right), \] entonces \(y\) es solución si cumple \[ y^{(n)}(t) = f\left(t, y(t), y'(t), \ldots, y^{(n-1)}(t)\right)\qquad \forall\, t\in I. \]

Nota: La definición de solución requiere que su dominio sea un intervalo \(I\). Esto tiene dos implicaciones fundamentales: 1. Conexión: No se trabaja en conjuntos no conexos. 2. Independencia de ramas: Una solución puede generar múltiples soluciones si su dominio no es conexo. Cada intervalo conexo del dominio define una solución diferente.

Soluciones explícitas e implícitas

Las soluciones se pueden expresar de forma explícita o implícita:

Definición 1.6 (Solución Explícita) Una solución es explícita si se puede escribir como una función de la variable independiente \(y = \phi(t)\).

Definición 1.7 (Solución Implícita) Una solución es implícita si se expresa mediante una relación \(G(t, y) = 0\) que define a una o varias funciones \(y=\phi(t)\) diferenciables en algún intervalo.

Ejemplo 1.7 La ecuación \(y' = -\dfrac{x}{y}\) tiene como solución implícita la ecuación \[ y^2 + x^2 = C \quad C \in \mathbb{R} \]

En forma explícita, la solución es \[ y = \pm \sqrt{C - x^2} \quad C \in \mathbb{R} \]

Nota: se observa que la solución implícita no es una función, cada solución implícita puede definir varias soluciones explícitas.

Código

data_fig_sol = {
  const pts = [];
  const R = 4;
  const n = 100;

  // 1. Datos Solución Implícita (Círculo completo ordenado)
  for (let i = 0; i <= n; i++) {
    const t = (i / n) * 2 * Math.PI;
    pts.push({
      x: R * Math.cos(t),
      y: R * Math.sin(t),
      title: "1. Solución Implícita",
      group: "impl"
    });
  }

  // 2. Datos Soluciones Explícitas (Dos ramas separadas)
  for (let i = 0; i <= n; i++) {
    const x = -R + (2 * R * i) / n;
    const y = Math.sqrt(Math.max(0, R*R - x*x));
    
    // Rama positiva
    pts.push({
      x, y, 
      title: "2. Soluciones Explícitas", 
      group: "pos"
    });
    
    // Rama negativa
    pts.push({
      x, y: -y, 
      title: "2. Soluciones Explícitas", 
      group: "neg"
    });
  }
  return pts;
}

Plot.plot({
  height: 200,
  width: 700,
  marginLeft: 40,
  grid: true,
  x: {domain: [-4.5, 4.5], label: "x"},
  y: {domain: [-4.5, 4.5], label: "y", ticks: 3},
  facet: {
    data: data_fig_sol,
    x: "title"
  },
  marks: [
    Plot.frame(),
    Plot.line(data_fig_sol, {
      x: "x", 
      y: "y", 
      stroke: d => d.group === "impl" ? "#17a2b8" : (d.group === "pos" ? "#27ae60" : "#e74c3c"),
      strokeWidth: 3,
      z: "group" 
    }),
    Plot.ruleX([0], {strokeOpacity: 0.3}),
    Plot.ruleY([0], {strokeOpacity: 0.3})
  ]
})

Comprobar que una función es solución

Para comprobar que una función es solución de una EDO, es suficiente con sustituirla en la ecuación y verificar que se cumple.

En los siguientes ejemplos se muestra soluciones de ecuaciones y su comprobación.

La ecuación \(y' = \dfrac{1}{x^2}\) tiene como familia de soluciones \(y = -\dfrac{1}{x} + C\) con \(C \in \mathbb{R}\).

Para comprobar que es solución, se sustituye en la ecuación:

\[ \frac{d}{dx} \left( -\frac{1}{x} + C \right) = \frac{1}{x^2} \]

Sin embargo esta función no está definida en \(x=0\). Esto quiere decir que la expresión análitica genera dos soluciones distintas, una definida en \((-\infty, 0)\) y otra definida en \((0, \infty)\). Se denominan soluciones maximales y son objeto de estudio en cursos más avanzados. Además, son solución también en cualquier subintervalo de sus dominios.

Ejemplo 1.8 Dada la ecuación anterior \(y' = -\dfrac{x}{y}\), se comprueba toda función diferenciable \(y\) que sarisfaga la ecuación implícita \(y^2 = -x^2 + C\) es solución.

Derivando ambos lados de la ecuación implícita respecto a \(x\):

\[ y^2 = -x^2 + C \]

\[ \frac{d}{dx}(y^2) = \frac{d}{dx}(-x^2 + C) \]

\[ 2y\frac{dy}{dx} = -2x \]

\[ y' = -\frac{2x}{2y} = -\frac{x}{y} \]

Por lo tanto, la solución implícita es solución de la ecuación.

Ejemplo 1.9 Dada la ecuación anterior \(y' = -\dfrac{x}{y}\), se comprueba que la función explícita \(y = \pm \sqrt{-x^2 + C}\) es solución.

Sustituyendo directamente:

\[ \frac{d}{dx} \left( \pm \sqrt{-x^2 + C} \right) = -\frac{x}{\pm \sqrt{-x^2 + C}} = -\frac{x}{y} \]

Por lo tanto, la solución explícita es solución de la ecuación. Además, es importante definir el intervalo en el que la función es solución, en este caso es solución en cualquier subintervalo de su dominio \(x \in (-\sqrt{C}, \sqrt{C})\). Nótese que la función está definida en el intervalo cerrado \([-\sqrt{C}, \sqrt{C}]\) pero no es solución en los extremos porque no es diferenciable en esos puntos. Este es un buen ejemplo de como la función puede ser solución en un intervalo pero no en su frontera.

De esta reflexión surgirá la necesidad de una teoría que intente extender soluciones a su intervalo maximal de existencia, caracterizar cuando este existe y como se comporta la solución cerca de la frontera.

La ecuación diferencial puede estar definida por partes. Veamos como se comporta:

Ejemplo 1.10 Sea la EDO \[ y' = \begin{cases} \cos x, & x \in (0, \pi), \\ 0, & x \notin (0, \pi). \end{cases} \]

Se puede comprobar que \(\sin x\) es solución en el intervalo \((0, \pi)\) pero no es solución en ningún intervalo que contenga puntos fuera de \((0, \pi)\). Habitualmente, se definen las soluciones a tramos intentando que sean continuas, por ejemplo:

\[ y(x) = \begin{cases} C_1, & x \leq 0, \\ \sin x + C, & x \in (0, \pi), \\ C_2, & x \geq \pi. \end{cases} \]

Esta familia de soluciones define tres soluciones de la ecuación diferencial definidas respectivamente en \((-\infty, 0), (0, \pi), (\pi, \infty)\). Ya que la derivada en \(\{0, \pi\}\) no está definida. Por ejemplo, el límite por la izquierda en \(0\) es \(0\) y por la derecha \(\cos 0 = 1\).

En función del modelo que represente nuestro problema, puede ser interesante definir la solución a tramos de forma que sea continua o incluso diferenciable. En este caso, se puede imponer solo la continuidad y se obtiene \(C_1 = C\) y \(C_2 = C + \sin \pi = C\). Por lo tanto, la familia de soluciones continuas es: \[ y(x) = \begin{cases} C, & x \leq 0, \\ \sin x + C, & x \in (0, \pi), \\ C, & x \geq \pi. \end{cases} \]

Por definición no son solución en los puntos \(\{0, \pi\}\) pero la diferenciabilidad en esos puntos puede no ser necesaria a la hora de estudiar el fénomeno que se está modelizando. En análisis avanzado, este tipo de funciones que cumplen la ecuación en casi todo su dominio se conocen como soluciones en sentido débil o soluciones de Carathéodory.

Código

viewof c1_piece = Inputs.range([-2, 2], {value: 0, step: 0.1, label: "C₁ (x ≤ 0)"})
viewof c_piece = Inputs.range([-2, 2], {value: 0, step: 0.1, label: "C (0 < x < π)"})
viewof c2_piece = Inputs.range([-2, 2], {value: 0, step: 0.1, label: "C₂ (x ≥ π)"})

Código

{ 
  const x_min = -2;
  const x_max = 5;
  const pi = Math.PI;
  
  // Generar datos para cada tramo
  const points = [];
  const n = 100;
  
  // Tramo 1: x <= 0
  for(let i = 0; i <= n; i++) {
    const x = x_min + (0 - x_min) * (i/n);
    points.push({x, y: c1_piece, type: "1. x ≤ 0"});
  }
  
  // Tramo 2: 0 < x < pi
  for(let i = 0; i <= n; i++) {
    const x = 0 + (pi - 0) * (i/n);
    points.push({x, y: Math.sin(x) + c_piece, type: "2. 0 < x < π"});
  }
  
  // Tramo 3: x >= pi
  for(let i = 0; i <= n; i++) {
    const x = pi + (x_max - pi) * (i/n);
    points.push({x, y: c2_piece, type: "3. x ≥ π"});
  }

  // Puntos frontera para destacar continuidad/discontinuidad
  const boundaryPoints = [
    {x: 0, y: c1_piece, type: "1. x ≤ 0", marker: true},
    {x: 0, y: Math.sin(0) + c_piece, type: "2. 0 < x < π", marker: true},
    {x: pi, y: Math.sin(pi) + c_piece, type: "2. 0 < x < π", marker: true},
    {x: pi, y: c2_piece, type: "3. x ≥ π", marker: true}
  ];

  return Plot.plot({
    height: 300,
    width: 600,
    marginLeft: 40,
    grid: true,
    x: {label: "x", domain: [x_min, x_max]},
    y: {label: "y", domain: [-3, 3]},
    color: {
      domain: ["1. x ≤ 0", "2. 0 < x < π", "3. x ≥ π"],
      range: ["#17a2b8", "#27ae60", "#e74c3c"]
    },
    marks: [
      Plot.ruleY([0], {strokeOpacity: 0.3}),
      Plot.ruleX([0, pi], {strokeOpacity: 0.3, strokeDasharray: "4"}),
      
      // La función a trazos
      Plot.line(points, {x: "x", y: "y", stroke: "type", strokeWidth: 3}),
      
      // Puntos en las fronteras
      Plot.dot(boundaryPoints, {
        x: "x", y: "y", 
        fill: d => d.type === "1. x ≤ 0" ? "#17a2b8" : (d.type === "2. 0 < x < π" ? "#27ae60" : "#e74c3c"), 
        r: 4
      }),
      
      // Texto indicando si es continua
      Plot.text(
        [{
          x: x_max, 
          y: 2.8, 
          label: (Math.abs(c1_piece - c_piece) < 0.05 && Math.abs(c2_piece - c_piece) < 0.05) 
            ? "CONTINUA" 
            : "DISCONTINUA"
        }], 
        {
          text: "label", 
          textAnchor: "end", 
          fill: (Math.abs(c1_piece - c_piece) < 0.05 && Math.abs(c2_piece - c_piece) < 0.05) ? "#27ae60" : "#e74c3c",
          fontWeight: "bold",
          fontSize: 14
        }
      )
    ]
  })
}

Modelos

Los modelos matemáticos que estudiaremos cumplen tres propiedades fundamentales:

Determinismo: el estado futuro y pasado del sistema está completamente determinado por su estado inicial. No existe ambigüedad: las soluciones determinadas por las condiciones iniciales son únicas.
Diferenciabilidad: las funciones que describen la evolución del sistema son continuas y derivables. Esta hipótesis, nos permitirá garantizar la existencia de soluciones.
Finito-dimensionalidad: el sistema se describe mediante un número finito de variables. Cada estado del sistema puede representarse como un punto en un espacio euclídeo de dimensión \(n\).

En conjunto, estas tres propiedades definen el marco clásico de las ecuaciones diferenciales ordinarias deterministas, que constituyen la base teórica de buena parte de la modelización matemática en ciencias e ingeniería.