Ecuaciones Diferenciales Ordinarias (EDOs)

Resolución simbólica y 
algunas interpretaciones gráficas

Realizado por el Prof. Pedro González Rodelas
Dpto. de Matemática Aplicada. Univ. de Granada

Fecha de la última revisión: 14/04/2015

Cualquier notificación de errata o sugerencia puede ser enviada al email: prodelas@ugr.es

Paquetes externos necesarios (usando la versión 3.0 de Mathematica)

<< Graphics`ImplicitPlot`

<< Graphics`PlotField`

Resolución simbólica de EDOs

En esta sección nos centraremos en la posible resolución analítica, aprovechando las capacidades de cálculo simbólico que proporciona Mathematica, de algunos tipos de Ecuaciones Diferenciales Ordinarias (EDOs), tanto lineales como no lineales. 
	Pero hay que dejar muy claro en este punto que, salvo las EDOs lineales y muy pocos tipos de las No Lineales, la mayoría de las Ecuaciones Diferenciales son bastante difíciles de resolver analíticamente, a no ser que se consiga identificarlas y clasificar como de alguno de los tipos que ya se tiene bien estudiados y se sabe como resolver, mediante técnicas específicas y algo elaboradas, en alguno de los casos. En todos los demás casos, casi siempre nos tendremos que conformar con realizar un pequeño estudio teórico previo relacionado con la posible existencia y unicidad de solución del problema en cuestión, seguido de una resolución numérica para calcular al menos una aproximación lo suficientemente buena del problema.

Mediante el comando DSolve

Mathematica dispone del comando DSolve para intentar la resolución simbólica de la Ecuación Diferencial que le indiquemos, ya sea sin o con condiciones iniciales o de contorno añadidas. Veámos a continuación cuál sería su sintaxis genérica.

?DSolve

Como vemos en la descripción general que nos indica la ayuda de Mathematica, este comando necesitaría tres tipos de argumentos bien diferenciados para su buen funcionamiento:  
	
	- en "eqn" deberíamos indicar tanto la(s) ecuacion(es) diferencial(es) que queramos resolver como la(s) condición inicial(es) y de contorno que hubiera en nuestro problema, todo ello entre llaves y separado por comas;
	
	- en "y" tendremos que indicar cuál es realmente nuestra(s) función(es) incógnita, dependiente(s) de la variable(s) independiente(s) correspondiente(s), con la letra y formato que aparezca en la(s) ecuación(es) diferencial(es) del argumento anterior; en caso de que sean varias habría que incluirlas entre llaves y separadas por comas, como suele ser habitual en Mathematica.
	
	- en "x" indicaríamos simplemente la(s) letra(s) que hayamos usado para la(s) variable(s) independiente(s) de nuestro problema, ya sean una o varias.

Veámos algunos ejemplos, tanto de  Ecuaciones diferenciales genéricas, como de Problemas de Valores Iniciales (PVIs) o de Contorno (PCs):

Ecuaciones Diferenciales en general

Resolvamos por ejemplo la EDO  y' - y = (e^x); nótese que dentro de la sentencia DSolve de hay que especificar la dependencia de la variable dependiente de la independiente y que la ecuación diferencial hay que indicarla con doble igual, como cualquier otra ecuación en  Mathematica. Nótese también que aunque hayamos escrito en el texto la función exponencial con la letra (e^x) , dentro de la orden de Mathematica no se nos puede olvidar usar la letra especial que hay dentro de la paleta de símbolos, o bien cambiar la letra "e" minúscula por la mayúscula, que ya si sería correctamente interpretada.

DSolve[y'[x] - y[x] == (E^x), y[x], x]

obteniendo la solución genérica y(x) = x (e^x) + (Subscript[c, 1]) (e^x).  Nótese  que Mathematica  irá denotando C[1], C[2], ... las distintas constantes de las que dependa la solución, dependiendo casi siempre del orden de la correspondiente ecuación diferencial de la que procede.

Veámos ahora otro ejemplo, con una ecuación diferencial de segundo orden     y'' = y + y' + sen x,    en la que la solución general va a depender de dos constantes en vez de sólo una

DSolve[y''[x] == y[x] + y'[x] + Sin[x], y[x], x] // Simplify

Conviene realizar una cierta simplificación (a través de la orden de Mathematica  Simplify, que podemos aplicar fácilmente tras la orden DSolve con la doble barra inclinada //).

Veámos ahora otro ejemplo de EDO no lineal  (nótese la expresión x(((t))^3) en uno de los términos de la ec. dif.) 
x'(t) + t x(t) = (t^3)x(((t))^3), en la que la solución general tiene varias ramas, además de depender de una constante (Subscript[c, 1]),   

		x(t) = ±(1/(Sqrt[1 + (t^2) - 2 (Subscript[c, 1]) (e^((t^2)))]))

DSolve[x'[t]+t x[t]==(t^3) x[t]^3,x[t],t]

Problemas de Valores Iniciales (PVIs)

En estos casos se usaría el comando DSolve  con el siguiente formato  (asumiendo que la variable dependiente es "y" y la independiente "x", aunque otras muchas posibilidades se podrían dar)

               DSolve[{ecuación diferencial, condición inicial}, y[x], x]
               
	Resolvamos por ejemplo el  PVI que se obtiene al añadir la siguiente condición incial y(0) = 1  a la ecuación diferencial resuelta anteriormente.

sol = 
      DSolve[{y'[x] - y[x] == (E^x),  y[0] == 1}, y[x], x]

Podemos ahora dibujar fácilmente la solución obtenida mediante el comandoPlot, indicándole un cierto intervalo de la variable independiente. Nótese que hemos aprovechado el comando de asignación que posee Mathematica  "/." y que permite asignarle a la función "y[x]" de la parte izquierda la solución obtenida anteriormente y que habíamos denominado "sol". El hecho de indicar "sol[[1]]" es porque hemos escogido la primera y única solución; si hubiera habido dos o más soluciones, éstas vendrían numeradas "sol[[2]]", "sol[[3]]", ...

Plot[y[x] /. sol[[1]], {x, 0, 2}]

Nótese que lo que hemos hecho es simplemente añadir dicha condición inicial en forma de nueva ecuación (con doble igual y usando corchetes en vez de paréntesis, ya que en Mathematica  así se indica cualquier argumento de una función) dentro de la orden DSolve.  Así hemos obtenido ya una solución concreta y única  y(x) = (e^x) (1 + x) , que ya no depende de ninguna constante. En realidad es cómo si hubiésemos calculado el valor que debe tomar la constante (Subscript[c, 1]) en el caso anterior para que se cumpla la condición inicial impuesta; es decir

	1= y(0) = 0 (e^0) + (Subscript[c, 1]) (e^0) = (Subscript[c, 1]).

Resolvamos ahora un PVI asociado a la EDO no lineal vista anteriormente, donde al imponer la condición inicial concreta ya se obtiene una única solución y se descarta la otra rama de posibles soluciones. Para realizar la  gráfica de la correspondiente solución usamos la misma orden de asignación que antes "/." pero si no nombramos a la solución podemos simplemente usar el resultado justo anterior, con el símbolo del porcentaje "%".

DSolve[{x'[t]+t x[t]==(t^3) x[t]^3,x[0] == (1/2)},x[t],t]Plot[x[t] /. %[[1]], {t, 0, 1}]

Problemas de Contorno (PCs)

La sintaxis general sería en este caso (asumiendo "y" y "x" como variables dependiente e independiente, respectivamente)
	
        DSolve[{ecuación diferencial, condiciones de contorno}, y[x],x]

Veámos ahora un ejemplo de resolución de un problema de contorno,  que al tener dos condiciones adicionales, tiene que venir a partir de una ecuación diferencial de segundo orden. 

			y'' = y + y' + sen x,  
      x ∈ [0, π], 
y(0) = 0,  
      y(π) = 1.

sol = 
      DSolve[{y''[x] == y[x] + y'[x] + Sin[x], y[0] == 0, 
            y[π] == 1}, y[x], x] // Simplify

Podemos ahora realizar la correspondiente gráfica, indicando con la opción "PlotRange->All" que se dibuje todo el rango posible de la variable dependiente en el intervalo indicado de la variable independiente.

Plot[y[x] /. sol[[1]], {x, 0, π}, PlotRange -> All]

EDOs de primer orden

Ecuación lineal de primer orden

Una ecuación diferencial ordinaria se dice lineal de primer orden, si se puede escribir de la siguiente forma

		y'(x) + a(x) y(x) = b(x)  
		
con a(x) y b(x) sendas funciones continuas, definidas en cierto intervalo I de la recta real.

En el caso particular de que se anule b(x) = 0, ∀x∈I⊂ℝ, se dice que dicha EDO es homogénea y se puede resolver escribiendo  y'(x)/y(x) = -a(x)  e integrando ambos miembros, respecto a la variable independiente x, obteniéndose   ln(y(x)) = -A(x), siendo A'(x) = a(x) una primitiva. Por tanto, cualquier función  y ≡ y(x)=C(e^(-A(x))) (con C∈ℝ, cte.) será una solución genérica de la ecuación homogénea. Y bastará con imponer una cierta condición adicional, del tipo y((Subscript[x, 0])) = (Subscript[y, 0])  por ejemplo, para determinar de forma unívoca una determinada solución, que satisfará el correspondiente PVI.

Resolvamos por ejemplo el siguiente PVI lineal homogéneo:  y' = -k y, y(0) = 1, que puede provenir de un cierto problema de transmisión de calor (según la ley de calentamiento/enfriamiento de Newton), un modelo Malthusiano de crecimiento de poblaciones, o bien de un típico problema de desintegración radioactiva.

DSolve[{y'[x] == -k y[x], y[0] == 1}, y[x], x]

En el caso no homogéneo, también resulta fácil comprobar que la solución de la EDO lineal de primer orden se puede expresar de la siguiente forma

		y ≡ y(x)=(c+B(x))(e^(-A(x))) (con c∈ℝ, cte.)

siendo A(x) cualquier primitiva de  a(x)  y  B(x) una primitiva de  b(x)(e^(A(x))).

\(
Clear[x, y, c, B]\)y[x_] := (c + B[x]) (E^(-A[x]))

y'[x] + a[x] y[x] - b[x]      %/.{A'[x]->a[x]}%/.{(E^(-A[x])) B'[x]->b[x]}

Resolvamos por ejemplo el siguiente PVI lineal no homogéneo:  y' +2x y = x, y(0) = 1.

Clear[x, a, b, A, B]a[x_] := 2  xb[x_] := x

A[x] = ∫a[x] \[DifferentialD]xB[x] = ∫b[x] (E^A)[x] \[DifferentialD]x

y[x] = (c + B[x]) (E^(-A[x]))y[x]

\(x0 = 0; \ y0 = 1; 
y[x] /. x -> x0\)Solve[% == y0, c]y[x] /. %[[1]]

Y comprobemos también el resultado mediante la orden DSolve

Clear[x, y]DSolve[{y'[x] + a[x] y[x] == b[x], y[x0] == y0}, y[x], x]

Ecuaciones no lineales de orden uno (caso no trivial)

Cuando se trata de resolver EDOs no lineales, ya la cuestión no resulta tan fácil, tanto para nosotros como para Mathematica. Intentemos resolver por ejemplo el siguiente PVI para la EDO no lineal de primer orden

	y'(x) = - (((((y(x))/x)) )^2), y(0) = 1
	
que podemos empezar a transformar de la siguiente forma  ((y' (x))/((y(x))^2)) = -(1/(x^2)), que integrando quedaría (1/(y(x)))=((-1)/x) + c, o bien y(x) = (x/(c x - 1)) .

	Pero esta familia de soluciones no cumpliría la condición inicial especificada, sino y(0)=0. De todas formas, ninguno de estos problemas se podría resolver sin complicaciones mediante la orden DSolve. Veámoslo a continuación.

DSolve[{y'[x] == -(((y[x]/x))^2), y[0] == 1}, y[x], x]

\(
Clear[x, y]\)solgen = 
      DSolve[{y'[x] == -(((y[x]/x))^2), y[0] == 0}, y[x], x]

Si ahora intentamos resolver la ecuación diferencial sin condición inicial alguna, Mathematica nos proporciona una solución genérica como la obtenida anteriormente de forma manual, que tampoco podríamos afirmar que es la única posible.

\(
Clear[x, y]\)solgen = DSolve[y'[x] == -(((y[x]/x))^2), y[x], x]

y[x] /. solgen[[1]]% /. {x -> 0}

Y vemos que podemos encontrar por ejemplo infinitas soluciones del PVI con la condición inicial y(0) = 0.

y[x_] := (x/(-1 + x C[1]))Plot[{y[x] /. C[1] -> 0, y[x] /. C[1] -> -1, 
        y[x] /. C[1] -> -2, y[x] /. C[1] -> -3}, {x, 0, 1}]

Respecto a la otra condición inicial que buscábamos, y(0) = 1, no hemos podido concluir nada hasta el momento. Vemos pues claramente las dificultades con las que nos podemos encontrar al tratar con ecuaciones diferenciales de tipo no lineal, necesitando en todo caso realizar un estudio previo acerca de la existencia y unicidad de solución, antes de poder afirmar nada con rotundidad.

Caso de variables separadas o separables

Consideremos la EDO  ((y^2) + x (y^2))y' + (x^2) - y (x^2) =  0, que de forma equivalente se puede expresar también  
(y^2)(1+x)y'   =  (x^2)(y - 1)  o incluso  ((y^2)/(y - 1))y' =(((x^2))/(1 + x)). 
	Este ejemplo nos muestra la idea subyacente detrás del tipo de ecuaciones diferenciales de variables separadas, ya que de alguna manera hemos conseguido expresar dicha ecuación mediante una igualdad en la que en cada uno de los miembros sólo aparece una de las variables. Para resolver ahora dicha ecuación diferencial, procederemos a integrar ambos miembros, cada uno respecto a la variable involucrada.

expry = ∫((y^2)/(y - 1)) \[DifferentialD]y // Simplify

exprx = ∫((x^2)/(1 + x)) \[DifferentialD]x // Simplify

Expand[exprx - expry] == Cte

Caso de ecuaciones homogéneas

Consideremos ahora la EDO  4x-3y+(2y-3x)y' = 0, que tiene la particularidad de poder ser expresada como  
y' = ((3 y-4  x)/(2 y - 3 x))=((3 (y/x)-4)/(2 (y/x) - 3)) y por ello puede considerarse como homogénea.

Bastará pues realizar el cambio de variable  u≡(y/x),  de manera que   
y' = u' x+ u = ((3 u-4)/(2 u - 3))  y podemos expresar la correspondiente ecuación diferencial para la nueva variable dependiente   u' ≡(du/dx)= (1/x)(((3 u-4)/(2 u - 3))-u)=(1/x)(((-2 (u^2)+6u-4)/(2 u - 3))), que se puede resolver separando variables e integrando como antes

			∫(1/x) dx=∫((2 u - 3)/(-2 (u^2)+6u-4))du+Cte

∫(1/x)\[DifferentialD]x-∫((2 u - 3)/(-2 (u^2)+6u-4))\[DifferentialD]u==Cte

Equivalentemente, 

	2  Log (x) +  Log (2 - 3 u + (u^2))-2Cte = Log(K (x^2)(2 - 3 u + (u^2))) = 0

o bien  \(TraditionalForm
      \`\(1 = \(K\ \(\(x\^2\)(2 - 3\ u + u\^2)\)\  = \ 
          \(K\ \(\(x\^2\)(2 - 3\ y\/x + \((y\/x)\)\^2)\)\  = \ 
            K\ \((2  x\^2 - 3\ x\ y + y\^2)\)\)\)\ \)\)

\(\(y\_+\)[x_] := 
      \(3  x\)\/2 + \@\(9  x\^2 - 4 \((2  x\^2 - 1\/K)\)\)\/2\)\(\(y\_-\)[x_] := 
      \(3  x\)\/2 - \@\(9  x\^2 - 4 \((2  x\^2 - 1\/K)\)\)\/2
\)\(K = 1; 
Plot[{\(y\_+\)[x], \(y\_-\)[x]}, {x, \(-10\), 10}]\)

Realicemos a continuación algunas comprobaciones

\(TraditionalForm
    \`K\ \((2  x\^2 - 3\ x\ \(y\_+\)[x] + \(y\_+\)[x]\^2)\)\  // Simplify\)

\(∂\_x\ \(y\_+\)[x]\)-((\(3 \(y\_+\)[x]\)-4  x)/(\(2 \(y\_+\)[x] - 3 x\)))% // Simplify

\(∂\_x\ \(y\_-\)[x]\)-((\(3 \(y\_-\)[x]\)-4  x)/(\(2 \(y\_-\)[x] - 3 x\)))% // Simplify

Caso de ecuaciones exactas

Consideremos la EDO   (2 x y - (x^2)) dx  + ((x^2) - 2  y) dy = 0.  Empezaremos comprobando si efectivamente se trata de una ecuación diferencial exacta

m[x_, y_] := 2  x y - (x^2)n[x_, y_] := (x^2) - 2  y(∂_y) m[x, y] == (∂_x) n[x, y]

Para obtener ahora la solución, procederemos de la manera habitual cuando se trata de este tipo de ecuaciones diferenciales

f[x, y] = ∫m[x, y] \[DifferentialD]x + h[y]

Solve[(∂_y) f[x, y] == n[x, y], h'[y]]h[y] = ∫(h'[y] /. %[[1]]) \[DifferentialD]y

f[x, y]

Así pues, la solución implícita de la ecuación diferencial sería 

			-((x^3)/3) + (x^2) y - (y^2) = K   (con  K∈ℝ, cte.)
y, podríamos definir como en algún caso anterior

\(\(y\_+\)[x_] := x\^2\/2 + \@\(x\^4 - 4 \((x\^3\/3 + K)\)\)\/2\)\(\(y\_-\)[x_] := x\^2\/2 - \@\(x\^4 - 4 \((x\^3\/3 + K)\)\)\/2
\)

NRoots[(x^4) - 4 (((x^3)/3) + K) == 0, x]

\(x1 = \(-1.16921063795860136`\); 
x2 = 1.90863201469387568`; \)

\(K = 1; 
graf1 = Plot[{\(y\_+\)[x], \(y\_-\)[x]}, {x, \(-10\), x1}]\)\(graf2 = Plot[{\(y\_+\)[x], \(y\_-\)[x]}, {x, x2, 10}]\)

Show[graf1, graf2]

Realicemos a continuación algunas comprobaciones

\(
\((\(-\(x\^3\/3\)\) + x\^2\ \(y\_+\)[x] - \(y\_+\)[x]\^2 - K)\) // 
      Simplify
\)\(\((\(-\(x\^3\/3\)\) + x\^2\ \(y\_-\)[x] - \(y\_-\)[x]\^2 - K)\) // 
      Simplify\)

\(∂\_x\ \(y\_+\)[x]\)-((2x \(\(y\_+\)[x]\)-(x^2))/(\(2\ \(y\_+\)[x]\)-(x^2)))% // Simplify

\(∂\_x\ \(y\_-\)[x]\)-(((x^2)-2x \(\(y\_-\)[x]\))/(\(x\^2 - 2\ \(y\_-\)[x]\)))% // Simplify

Uso de factores integrantes

Consideremos la EDO  (y - 8x (y^2))dx+(3x-16(x^2) y)dy = 0, que podría ser o no exacta. Lo comprobaremos y en caso negativo, la técnica del factor integrante trata de proporcionarnos un factor multiplicativo que ayude a hacerla de ese tipo.

m[x_, y_] := y - 8  x (y^2)n[x_, y_] := 3  x - 16  (x^2) y(∂_y) m[x, y] == (∂_x) n[x, y]

Vemos claramente que no lo es, pero vamos a comprobar si admitiría un factor integrante de la forma μ ≡ μ(x) , para lo cual la siguiente expresión debería depender tan sólo de la variable independiente x.

exprx = 
      Simplify[(((∂_y) m[x, y] - (∂_x) n[x, y])/n)[x, y]]

Probemos ahora con un factor integrante de la forma μ ≡ μ(y) , para lo cual la siguiente expresión debería depender tan sólo de la variable independiente y.

expry = 
      Simplify[(((∂_y) m[x, y] - (∂_x) n[x, y])/m)[x, y]]

En este caso, sí que ha habido suerte y podemos buscar el factor integrante μ ≡ μ(y)   gracias a la fórmula siguiente

\(\(
μ[y] = Exp[\(-\(∫expry \[DifferentialD]y\)\)]\)\)

μ[y]

my[x_, y_] := m[x, y] μ[y]ny[x_, y_] := n[x, y] μ[y](∂_y) my[x, y] - (∂_x) ny[x, y]% // Simplify

Para obtener ahora la solución de la ecuación diferencial original, empezaremos calculando

expry = ∫my[x, y] \[DifferentialD]x

y definimos la función

Clear[x, y, f, h]f[x_, y_] := expry + h[y]

f[x, y]

y resolvemos finalmente la ecuación

Solve[(∂_y) f[x, y] == ny[x, y], h'[y]]h[y] = ∫(h'[y] /. %[[1]]) \[DifferentialD]y

f[x, y]

Siendo finalmente la solución de la ecuación diferencial la dada implícitamente mediante la expresión   x (y^3) - 4  (x^2) (y^4) =  k, para k∈ℝ  (cte).

Otras ecuaciones de primer orden (Ec. de Bernouilli)

La ecuación diferencial ordinaria  y'(x) + P(x) y(x) = Q(x) ((y(x))^n)  se conoce como ecuación de Bernouilli. Aunque se puede ver con claridad que no es lineal, sí que puede reducirse a una de este tipo con el siguiente cambio de variable dependiente  u(x)≡ ((y(x))^(1 - n)), transformándose en 

	u'(x) =(1-n) ((y(x))^(-n)) y'(x) = (1-n) ((y(x))^(-n)) (Q(x) ((y(x))^n) - P(x) y(x)) = (1-n) (Q(x)  - P(x) u(x))

Resolvamos el siguiente ejemplo:    y'(x) + ((y(x))/x) = (x^3) ((y(x))^3)

Clear[x, P, Q]P[x_] := (1/x)Q[x_] := (x^3)

\(n = 3; 
ecdifBernouilli = \(y'\)[x] + P[x] y[x] == Q[x]\ y[x]\^n\)

ecdifmod = u'[x] == (1 - n) (Q[x] - P[x] u[x])

solEDu = DSolve[ecdifmod, u[x], x]

solEDu /. {u[x] -> y[x]^(1 - n)}

\(\(y\_+\)[x_] := \((\(-x\^4\) + x\^2\ C[1])\)\^\(-2\)\)\(\(y\_-\)[x_] := \(-\((\(-x\^4\) + x\^2\ C[1])\)\^\(-2\)\)\)

\(ϵ = 0.1; \ c = 1; 

    Plot[Evaluate[{\(y\_-\)[x], \(y\_+\)[x]} /. C[1] -> c], {x, ϵ, 
        \@c}]\)

EDOs de segundo orden

Lo primero que hay que aclarar, lo mismo que ya hicimos en el caso de ecuaciones de primer orden no lineales, es que no existen métodos generales de resolución analítica de EDOs de segundo orden, salvo algunos casos muy concretos que iremos indicando.

Caso de coeficientes constantes

Se trata de ecuaciones diferenciales de la siguiente forma

		y''(x) + a y'(x) + b y(x) = g(x)
		
distinguiendose dos casos bien diferenciados, como en el caso de ecuaciones de primer orden:
               - el homogéneo (cuando g(x)≡ 0) y 
               - el no homogéneo (en caso contrario).

En el caso homogéneo, supondremos soluciones de la forma  y(x) ≡  (e^(m x)) y saquemos factor común al introducir ααlas derivadas en la correspondiente ecuación diferencial, obteniéndose que
(e^(m x))((m^2) + a m + b)= 0   para todo x de la variable independiente donde se tenga definida la ec. dif.
Por ello, dado que la función exponencial no se anula, acabaremos teniendo la siguiente ecuación polinómica  (m^2) + a m + b = 0 para el coeficiente del exponente. Esta es la denominada ecuación auxiliar, polinómica de segundo orden, y según los valores del correspondiente discriminante  d = (a^2) - 4  b podemos distinguir varias posibilidades:

a) discriminante positivo: en cuyo caso se tienen dos soluciones reales distintas de la ec. auxiliar,  m=α ≡((-a)/2)+((Sqrt[d])/2) y  m=β≡((-a)/2)-((Sqrt[d])/2), siendo la solución general de nuestra ecuación diferencial no homogénea

		y(x) = \(TraditionalForm
      \`c\_1\ e\^\(α\ x\) + c\_2\ e\^\(β\ x\)\) +\(TraditionalForm
      \`\(1\/\(α - β\)\) 
        \((\(e\^\(α\ x\)\) \(P(x)\) - \(e\^\(β\ x\)\) \(Q(x)\))
          \)\), con (Subscript[c, 1]), (Subscript[c, 2])∈ℝ constantes.
para P(x) una primitiva de la función (e^(-α x)) g(x)  y  Q(x)  una primitiva de  (e^(-β x)) g(x), siendo claramente ambos nulos si estuvieramos en el caso homogéneo.
		
b) discriminante negativo: en este caso se tienen dos soluciones imaginarias conjugadas de la ec. auxiliar,  m=α ≡((-a)/2)+i((Sqrt[-d])/2) y  m=β≡((-a)/2)-i((Sqrt[-d])/2), siendo la solución general de nuestra ecuación diferencial no homogénea (siendo P(x) y Q(x)  nulos en el caso homogéneo).

		y(x) = (e^(((-a)/2) x))((Subscript[c, 1]) cos(((Sqrt[-d])/2) x)+(Subscript[c, 2]) sen(((Sqrt[-d])/2) x))+(1/(α - β)) 
              ((e^(α x)) P(x) - 
                  (e^(β x)) Q(x)) , con (Subscript[c, 1]), (Subscript[c, 2])∈ℝ constantes.
		
c) discriminante nulo: teniéndose en este caso una única solución real doble de la ec. auxiliar,  m=α ≡((-a)/2), siendo la solución general de nuestra ecuación diferencial no homogénea

		y(x) =  
      (e^(α x))((Subscript[c, 1]) + (Subscript[c, 2]) x) + (e^(α x)) P(x) 
       , con (Subscript[c, 1]), (Subscript[c, 2])∈ℝ constantes.

Para la resolución ahora del caso homogéneo, sólo tendríamos que recordar que entonces P(x)  = 0,

Otras EDOs interesantes

Ecuación de Ricatti

Es una ecuación general del tipo    y'+(y^2) + P(x) y + Q(x) = 0
que no puede resolverse fácilmente, salvo en casos especiales. Su importancia radica en la relación que tiene con la ecuación lineal de segundo orden   z''+P(x) z' +Q(x) z =0, sin más que realizar el siguiente cambio y≡\(TraditionalForm\`z'\/z\). Por tanto, si y ≡ y(x) es solución de la ecuación de Ricatti, entonces z(x) ≡ C(e^(A(x) )) será solución de la correspondiente ecuación lineal de segundo orden, siendo A(x) una primitiva de y(x). Así, en caso de conocerse alguna solución particular (Subscript[y, p])(x) , entonces la solución general puede determinarse de la siguiente forma:
- escribimos  y(x) ≡ (Subscript[y, p])(x) + (1/(v(x)))
- derivamos ahora y sustituimos, quedando que  v' - (P(x)+2(Subscript[y, p])(x))v(x) = 1
- si conseguimos ahora R(x) una primitiva de P(x)+2(Subscript[y, p])(x), bastará con resolver ahora la ecuación resultante, obteniendose que v(x) = (e^(R(x)))(S(x) + A(x)), para S(x) una primitiva de (e^(-R(x))).

Ecuación de Euler lineal

Es una ecuación lineal de orden n que se puede escribir como 

	\(TraditionalForm
      \`\(x\^n\ y\^\(n)\) + \(a\_1\) x\^\(n - 1\)\ y\^\(n - 1)\) + ...\) + 
          \(a\_\(n - 1\)\) x\ y' + a\_n\ y = g(x)\)
	
para (Subscript[a, 1]),...,(Subscript[a, n - 1]),(Subscript[a, n]) constantes reales. Y que puede reducirse fácilmente a otra de tipo lineal, pero ahora con coeficientes constantes, sin más que realizar el cambio de variable independiente x ≡ (e^t), pudiéndose escribir la nueva ecuación de la siguiente forma

	\(TraditionalForm
      \`\(y\^\(n)\) + \(q\_1\) y\^\(n - 1)\) + ...\) + q\_\(n - 1\)\ y' + 
          q\_n\ y = g(e\^t)\)

con  (Subscript[q, 1]),...,(Subscript[q, n - 1]),(Subscript[q, n]) constantes reales también.

Veamos un ejemplo concreto, resolvamos la ecuación diferencial  (x^2) y'' + 6  x y' + 6  y = (1/(x^2)).

ecdifEuler = 
      (x^2) y''[x] + 6  x y'[x] + 6  y[x] == (1/(x^2))

sol = DSolve[ecdifEuler, y[x], x]

ecdifEulerbis=ecdifEuler//.{x -> (E^t),y[(E^t)] -> w[t],(E^t) y'[(E^t)]->w'[t],(E^(2 t)) y''[(E^t)]->w''[t]}

solbis = DSolve[ecdifEulerbis, w[t], t]

Simplify[w[t] /. solbis[[1]]]% /. t -> Log[x]

Algunas Interpretaciones geométricas

Si tenemos que resolver, o al menos estudiar, la ecuación diferencial    y'(x) = f(x, y(x))  y sabemos que tendrá solución única para cada condición inicial y(a)  = α, podríamos representar gráficamente cada una de las curvas solución, para cada valor del parámetro α.

Veamos un ejemplo concreto,  y'(x) = x - y(x), cuya solución general se expresa
y(x) = -1 + x + (Subscript[c, 1]) (e^(-x)), como podemos comprobar fácilmente usando el comando DSolve.

solgen = DSolve[y'[x] == x - y[x], y[x], x]

Mediante la orden Evaluate podemos realizar una representación gráfica de una tabla con las soluciones obtenidas para cada uno de los valores de la constante en un cierto rango.

grafsolgen = 
      Plot[Table[
            (y[x] /. solgen[[1]]) /. C[1] -> c, {c, -1, 1, 0.2}]
             // Evaluate, {x, 0, 5}, PlotStyle -> RGBColor[1, 0, 0]]

Campo de direcciones

En general, para cada par de valores (a,α)  podemos saber de antemano cuál será la pendiente de la única curva solución  (Subscript[y, α])(x)  de la ecuación diferencial  que pasa por dicho punto, sin más que evaluar la función de dos variables que define dicha ecuación; es decir  y'(a) = f(a, α)  será la correspondiente pendiente. Una técnica habitual en el estudio del comportamiento de las soluciones de una ecuación diferencial en cierta región consiste en representar en cada punto del plano un pequeño vector cuya pendiente coincida con la pendiente así obtenida. A esta representación gráfica se la denomina campo de direcciones de la ecuación diferencial.

Mathematica dispone de la orden PlotVectorField, incluida dentro del paquete Graphics`PlotField`  para representar el correspondiente campo de direcciones.

?PlotVectorField

Empecemos definiendo la función de dos variables que define la ecuación diferencial correspondiente

Clear[x, y, f]f[x, y] := x - y

?f

\(\(x0 = 0; \ x1 = 5; \ y0 = \(-2\); \ y1 = 4; 

    campodir = 
      PlotVectorField[{1, x - y}, \ {x, \ x0, \ x1}, \ {y, \ y0, \ y1}]\ 
    \)\)

Podemos ahora representar conjuntamente este campo de direcciones y la familia de soluciones obtenidas anteriormente, para distintos valores de la condición inicial y vemos claramente cómo dicho campo de direcciones nos indica de una manera gráfica como se van a comportar las correspondientes soluciones de la ecuación diferencial que pasan por cada punto del plano.

Show[grafsolgen, campodir]

Un segundo ejemplo

Clear[x, y, f]f[x_, y_] := Sin[x + y]

?f

solgen = DSolve[y'[x] == f[x, y[x]], y[x], x]

Mathematica nos avisa de la posibilidad de que haya más soluciones, aparte de las obtenidas mediante el comando DSolve. De hecho puede probarse por otros medios que por cada punto del plano pasa una y sólo una solución de la ecuación diferencial, y vamos a representar el correspondiente campo de direcciones en un determinado rectángulo para hacernos una idea de como pueden ser dichas soluciones.

\(\(x0 = \(-2\) π; \ x1 = 2  π; \ y0 = \(-2\) π; \ 
    y1 = 2  π; 

    campodir = 
      PlotVectorField[{1, f[x, y]}, \ {x, \ x0, \ x1}, \ {y, \ y0, \ y1}, 
        PlotPoints -> 25]\ \)\)

Como sabemos, el campo de direcciones de una ecuación diferencial nos puede proporcionar una valiosísima información gráfica acerca de como pueden comportarse las posibles soluciones. En este caso concreto que hemos representado, cabe pensar que haya soluciones que sean rectas (siguiendo el campo de direcciones rectilíneo que se puede entrever) y que el resto de soluciones oscilan de entre varias de estas rectas.

Esta ecuación diferencial se puede convertir en una de variables separables, sin más que efectuar el siguiente cambio de variable  u≡ x+y, con lo que u' ≡ 1+y' =1+ sen(u), que podemos integrar ahora sin demasiado esfuerzo, ya que podemos escribir (du/(sen(u) + 1)) =  dx y bastará con calcular las primitivas correspondientes

(∫(1/(Sin[u] + 1)) \[DifferentialD]u // Simplify) == 
      ∫1 \[DifferentialD]x + k\(% /. {u -> x + y}
\)

\(ecsols = 
      Table[1 - Cos[x + y] + Sin[x + y] == \((x + k)\) \((1 + Sin[x + y])\), {
          k, \(-5\), 5, 0.5}]; 

    grafsols = 
      ImplicitPlot[ecsols, {x, x0, x1}, {y, y0, y1}, 
        PlotStyle -> RGBColor[1, 0, 0]]\)

Mostramos a continuación de manera conjunta estas gráficas y las del campo de direcciones asociado, para poder apreciar más claramente la influencia de dicho campo en las soluciones de la ecuación diferencial.

Show[campodir, grafsols]

Curvas isoclinas

Para una ecuación diferencial de primer orden  y'(x) = f(x, y(x)),  cada curva de nivel  f(x,y) = c (con c, constante) define lo que también se denomina "curva isoclina", de forma que se puede asegurar que en los puntos de una tal curva isoclina la pendiente de la solución de nuestra ecuación diferencial es siempre la misma (y'(x) =c, constante)

Vamos a considerar de nuevo la ecuación diferencial estudiada más arriba   y'(x) = sen(x+y(x))  y  implícitamente alguna de las curvas isoclinas asociadas.

\(curvasisoclinas = Table[f[x, y] == c, {c, \(-2\), 2, 0.2}]; 

    grafisoclinas = 
      ImplicitPlot[curvasisoclinas, {x, x0, x1}, {y, y0, y1}, 
        PlotStyle -> RGBColor[0, 1, 0]]\)

Representemos ahora conjuntamente, tanto estas curvas isoclinas como el correspondiente campo de direcciones y las soluciones obtenidas.

Show[campodir, grafsols, grafisoclinas]

Ecuación Diferencial de una familia de curvas

Si lo que nos planteamos ahora es el problema contrario al que hemos venido considerando hasta ahora: encontrar más bien una ecuación diferencial cuyas soluciones sean las curvas de una cierta familia n-paramétrica dada, tendremos entonces que derivar respecto de la variable independiente tantas veces como constantes haya y despejar dichas constantes, para obtener la ecuación diferencial o sistema de ecuaciones resultante.

Ejemplo

Supongamos que necesitamos hallar la EDO cuyas soluciones son las circunferencias cuyos centros están sobre el eje Ox  y son tangentes a las bisectrices de los cuadrantes primero y cuarto, cumpliéndose pues que   	
								(((x - (Subscript[x, c])))^2)+(y^2)=(r^(2 )) =  (((((Subscript[x, c])))^2)/2)  
ya que  r=dist(((Subscript[x, c]),0),{(x,y)∈(ℝ^2)  /    y - x = 0}) = \(TraditionalForm
      \`\(|\( x\_c | \)\)\/\@2\). 
	Así pues, derivando implícitamente la ecuación anterior respecto de la variable x, obtendríamos que  2(x - (Subscript[x, c]))  + 2  y(x) y'(x) =0,  que podríamos simplificar como simplemente x + y y' = (Subscript[x, c]) ; quedándo finalmente la siguiente ecuación diferencial, al sustituir esta expresión en la ecuación de la circunferencia del comienzo:
	
	  (x -((( x  +  y y' ))^2)+(y^2)= (((((x  +  y y')))^2)/2)=(1/2)((x^2)+2xyy'+(((yy'))^2))
	  
que simplificando ahora convenientemente, obtendríamos la siguiente ecuación diferencial asociada a las circunferencias consideradas:

			   (((y y' ))^2) + 2 (y^2) - 2xyy' = (x^2)

?Circle

\(x\_c = 5; θ\_min = π\/4; θ\_max = \(7  π\)\/4; 

    arco = Graphics[
        Circle[{x\_c, 0}, Abs[x\_c]\/\@2, {θ\_min, θ\_max}]]; 

    \(y\_+\)[x_] := x; \ \(y\_-\)[x_] := \(-x\)\)\(Plot[{\(y\_-\)[x], \(y\_+\)[x]}, {x, \(-10\), 10}]\)\(Show[%, arco, AspectRatio -> Automatic]
\)

Trayectorias ortogonales e isogonales

Los problemas relativos a trayectorias ortogonales no sólo son muy importantes en Matemática pura y aplicada, sino que también aparecen de forma natural en el estudio de muchos problemas de Física e Ingeniería. Por ejemplo, en la teoría de flujo de fluidos, a las familias ortogonales de curvas se las denominará líneas equipotenciales y líneas de corriente, respectivamente. En teoría del calor, éstas se denominarán líneas isotermas y líneas de flujo, etc.

Dada la ecuación de una familia de curvas uniparamétricas  F(x,y,C) = 0, diremos que una curva γ  es una trayectoria ortogonal, si ésta forma un ángulo recto con cada una de las curvas de la familia. Pero ¿cómo encontrar dichas curvas ortogonales, dada la familia de curvas? Simplemente se puede afirmar que  F(x,y,((-1)/y)') = 0 será la ecuación diferencial de dichas trayectorias ortogonales, ya que si   (Overscript[t, →]) es el vector tangente a una de las curvas de la familia y (Overscript[n, →])  resultase ser el correspondiente vector normal (girado 90∘  a partir del anterior), y por tanto el vector tangente a la curva ortogonal  γ ,  entonces la derivada correspondiente a esta trayectoria ortogonal se puede calcular a partir de la original, resultando ser ((-1)/y)'. El concepto de trayectorias isogonales viene a generalizar el concepto de curvas ortogonales, cuando se permite que el ángulo que forma dicha trayectoria con todas las curvas de la familia sea un cierto ángulo constante, pero no necesariamente de 90∘.

Así pues, si de antemano se sabe que la ecuación diferencial satisfecha por las soluciones de nuestra ecuación diferencial original es  y'(x) = f(x,y(x)), entoces bastaría con resolver la ecuación diferencial y'(x) = ((-1)/(f(x, y(x)))), para obtener las respectivas pendientes de las trayectorias ortogonales asociadas.

Ejemplo

Para hallar las trayectorias ortogonales a la familia  (x^2) + (((y - a))^2) = (a^2)
empezaremos derivando esta expresión respecto a la variable independiente x, para obtener la correspondiente Ecuación Diferencial de la familia de curvas, obteniéndose

	2x+2(y-a) y' = 0   con lo que a = y +(x/y)'  y sustituyendo esta expresión en la ecuación implícita de la familia original, nos queda  (x^2) + (((-(x/y)'))^2)=(((y +(x/y)'))^2)
y finalmente, simplificando adecuadamente:

		y' = 2((x y)/((x^2) - (y^2)))
		
y la ecuación diferencial que definirá las trayectorias ortogonales correspondientes será 

		y' = (((y^2) - (x^2))/(2 x y)).

Intentemos ahora realizar este mismo proceso de forma mecánica, con la ayuda del cálculo simbólico del que dispone Mathematica

Clear[x, a]eccurvas = (x^2) + (((y[x] - a))^2) == (a^2)

(∂_x) eccurvassol = Solve[%, a] // Simplify

eccurvas /. sol[[1]]Solve[%, y'[x]]

f[x_, y_] := ((2  x y)/((x^2) - (y^2)))

((-1)/f)[x, y]

solcurvas = DSolve[y'[x] == f[x, y[x]], y[x], x]

\(grafcurvas = 

      Plot[Evaluate[
			
          Table[
				{y[x] /. solcurvas[\([1]\)], 
                y[x] /. solcurvas[\([2]\)]} /. C[1] -> c, {c, 0.1, 1, 0.1}]

          				], {x, 0, 1\/2},  (*AspectRatio -> Automatic, *) 
			
        PlotStyle -> RGBColor[0, 0, 1]]\)

\(Table[{y[x] /. solcurvas[\([1]\)], y[x] /. solcurvas[\([2]\)]} /. 
        C[1] -> c
, {c, 1, 2}]\)

\(grafcurvas = 

      Plot[Evaluate[

          			{y[x] /. solcurvas[\([1]\)], y[x] /. solcurvas[\([2]\)]} /. 
            C[1] -> 1
			]
	, {x, 0, 1\/\(2  C[1]\) /. C[1] -> 1}, 
        AspectRatio -> Automatic, PlotStyle -> RGBColor[0, 0, 1]]\)

soltrayort = DSolve[y'[x] == ((-1)/f)[x, y[x]], y[x], x]

{y[x] /. soltrayort[[1]], y[x] /. soltrayort[[2]]} /. 
      C[1] -> c

\(graftrayort = 

      Plot[Evaluate[

          			{y[x] /. soltrayort[\([1]\)], y[x] /. soltrayort[\([2]\)]} /. 
            C[1] -> 1
				], {x, 0, 1}, AspectRatio -> 1]\)

\(graftrayort = 

      Plot[Evaluate[
			
          Table[
			{y[x] /. soltrayort[\([1]\)], 
                y[x] /. soltrayort[\([2]\)]} /. C[1] -> c, {c, 0.1, 1, 0.1}]

          				], {x, 0, 1}]\)

\(graftrayort = 

      Plot[Evaluate[
			
          Table[
			{y[x] /. soltrayort[\([1]\)], 
                y[x] /. soltrayort[\([2]\)]} /. C[1] -> c, {c, 0.1, 1, 0.1}]

          				], {x, 0, 1}, AspectRatio -> Automatic]\)

Show[grafcurvas, graftrayort, AspectRatio -> Automatic]

Curvas isotermas y curvas de flujo de calor

Si ahora la función T(x,y) representara la temperatura en el punto del plano (x,y), entonces la ecuación de las curvas de nivel T(x,y) = c (para c∈ℝ, una constante), se llamarán isotermas y las correspondientes trayectorias ortogonales serán las curvas de flujo de calor. Es posible que las que conozcamos de antemano sean precisamente estas últimas y necesitemos pues calcular a partir de ellas las correspondientes isotermas. En este caso, como sabemos que ambas son ortogonales entre sí, lo que tendremos que hacer es calcular las trayectorias ortogonales de la familia de curvas que conozcamos, en este caso las de flujo de calor.

Ilustremos el proceso con un ejemplo, en el que se supone que dicha familia de curvas para el flujo de calor viene dada por la ecuación  (y^2) + 2 x y - (x^2) = k  (para k∈ℝ, constante). Realizaremos todo el proceso con la ayuda de la capacidad de cálculo simbólico de Mathematica.

Clear[x, y, k]curvasflujo = y[x]^2 + 2  x y[x] - (x^2) == k

A continuación deberemos de calcular la ecuación diferencial asociada a esta familia de curvas, para lo cuál hay que derivar y eliminar la constante que aparece.

ecdif = D[curvasflujo, x]

Por último, para calcular la ecuación de las trayectorias ortogonales, que será la que nos proporcione las buscadas isotermas, lo que haremos será sustituir en esta ecuación diferencial  y'  por ((-1)/y)', ya que sabemos que el producto de las pendientes dos curvas ortogonales del plano, debe ser -1.

ecdifisotermas = 
      ecdif /. {y'[x] -> ((-1)/(y'))[x]} // Simplify

Resolviendo esta nueva ecuación diferencial obtendremos las curvas isotermas buscadas.

solisotermas = DSolve[ecdifisotermas, y[x], x]

Representemos gráficamente algunas de estas curvas isotermas, sin más que darle distintos valores a la constante de la que dependen.

\(curvasisotermas = 
      Table[{y[x] /. solisotermas[\([1]\)], y[x] /. solisotermas[\([2]\)]} /. 
          C[1] -> c, {c, 0.1, 1, 0.1}]; 

    grafisotermas = 
      Plot[Evaluate[curvasisotermas], {x, \(-5\), 5}, 
        PlotRange -> {\(-5\), 5}, AspectRatio -> Automatic]\)

\(curvasflujo = Table[y\^2 + 2  x\ y - x\^2 == k, {k, \(-10\), 10, 1}]; 

    grafsflujo = 
      ImplicitPlot[curvasflujo, {x, \(-5\), 5}, {y, \(-5\), 5}, 
        PlotStyle -> RGBColor[0, 0, 1]]\)

Show[grafisotermas, grafsflujo, AspectRatio -> Automatic]

Ejercicios

1.- Resolver, siempre que sea posible, las siguientes EDOs de primer orden

a)  2(x^2)+2(y^2)+(4x y+3(y^2))y'=0

b) 3(((y - 1))^2)y' = 2 + sen(x)

c) ((y^4) - (y^2) + 1) y' = (x^3) y  -  y

d) (2 + y) y' = 2x  + 3 y

e) y' = ((x  +  y)/(x - y))

f) y' = \(TraditionalForm
      \`\(2  x\ \(cos(y)\)\  + 3  x\^2\ y\)\/\(\(-x\^3\) + 
          \(x\^2\) \(sen(y)\) + y\)\)

g) (((y' + y) )^2) - (((x  -  y))^2) = 0

h) y = y' x + g(y')    (Ecuación de Clairaut)

i) y = y' P(x) + Q(x)(((y'))^n)

2.- Resolver, siempre que sea posible, las siguientes EDOs de segundo orden

a)  y''-2y'-3y=0

b)  y''+4y'+13y=0

c)  y''-4y'+4y=0

d)  y''-3y'+2y = cos(3x)

e)  y''+4y'+4y = cos(x)

f)  y'' + y= 0

3.- Resolver el siguiente PVI

a)  y''-y'-2y = 4(x^2),  y(0) = 0, y'(0) = 1.

4.- Representar el campo de direcciones asociado a las siguientes EDOs y estudiar su comportamiento y posibles soluciones a partir de este. Intentar también su resolución mediante la orden DSolve.

a)  y' = y\(TraditionalForm\`\^2\)+(e^(x + y)), para -5<x<5.

b)  y' = sen(x) -(y/x), para -5<x<5 (atención con los valores próximos a 0).

5.- Calcular y representar las trayectorias ortogonales de las siguientes familias de curvas.

a)  Las parábolas x = k (y^2),  con k∈ℝ  cte.

b)  Las parábolas y = k (x^2),  con k∈ℝ  cte.