Saltar al contenido

¿Cómo resolver la ecuación diferencial con la integral de Duhamel?

Después de tanto batallar hemos hallado la contestación de este enigma que muchos lectores de nuestro sitio web han tenido. Si deseas compartir algún detalle no dudes en dejar tu información.

Solución:

NDSolve actualmente no puede manejar este tipo de ecuación diferencial, LaplaceTransform es tu amigo. Dado que en este caso la transformada de Laplace inversa no se puede hacer analíticamente por InverseLaplaceTransform, necesita la ayuda del paquete de inversión numérico de Laplace además:

eq = 0.01 - 6.25 x[t] + (1.2 Integrate[x'[t - τ]/Sqrt[τ], τ, 0, t])/10^7 == 16 x''[t];
ic = x[0] == 0, x'[0] == 0;

f[s_] = With[l = LaplaceTransform, 
   Solve[l[Rationalize[eq, 0], t, s] /. Rule @@@ ic, l[x[t], t, s]]][[1, 1, -1]]

sol = FunctionInterpolation[GWR[f, t], t, $MachineEpsilon, 10]

Plot[sol[t], t, 0, 10]

ingrese la descripción de la imagen aquí

Observación:

  1. eq es Rationalized porque el primer argumento de GWR no debe contener números aproximados.

  2. $MachineEpsilon se elige como el límite inferior de FunctionInterpolation porque GWR tiene alguna dificultad para calcular la inversión en $ t = 0 $.

La figura de sol ya se ve bien, pero si no está satisfecho con la precisión actual de FunctionInterpolation, puede mejorarlo agregando derivadas de la expresión a su primer argumento. Esto se menciona en el Alcance del documento de FunctionInterpolation (¿Por qué no se coloca en un lugar más visible?):

gwr[t_?NumericQ] := GWR[f, t]

solOrder3 = 
 FunctionInterpolation[
  D[[email protected]
, t, #] & /@ Range[0, 3] // Evaluate, t, $MachineEpsilon, 10] (* Notice Chop is necessary in this case *) Plot[solOrder3[t] // Chop, t, 0, 10] (* Error check *) Plot[Subtract @@@ eq /. x -> solOrder3 // Abs // Evaluate, t, ##, PlotPoints -> 20, MaxRecursion -> 2, PlotRange -> All] & @@@ 0., 0.002, 0.002, 1, 1, 10 // GraphicsRow

ingrese la descripción de la imagen aquí

Publicación relacionada:

Cómo trazar y resolver la solución numérica de una ecuación integro-diferencial

Se puede usar la iteración de tipo Picard para obtener la solución: usando una aproximación para x'[t] (en la integral), podemos integrar la EDO para obtener una nueva aproximación. Sorprendentemente, converge en solo dos pasos. Mi pensamiento original fue avanzar paso a paso por la integración utilizando las herramientas de tutorial/NDSolveStateData para construir una interpolación de x'[t] en cada paso para su uso en el término de integración; eso resultó demasiado difícil de administrar (o tal vez lo había configurado de una manera que lo hacía difícil).

La aproximación de x'[t] está representado por xp[t]. Comenzamos con la suposición inicial para que sea xp[t] == 0.01 t, que corresponde a extrapolar de las condiciones iniciales (por inspección – uno podría resolver la EDO para x''). (En realidad, comenzando con xp[t] == 0 funciona casi tan bien y hace que la primera iteración sea más rápida). Ponemos el factor de integración en una función de caja negra separada y0. Sumando la ecuación algebraica ficticia y[t] == y0[t] al sistema ayuda con la precisión.

ClearAll[xp, y0, t, x, y];

xp = 0.01 # &;
y0[t_?NumericQ] := NIntegrate[xp[t - τ]/Sqrt[τ], τ, 0, t];

ode = 0.01 - 6.25 x[t] + 1.2 y0[t] / 10^7 == 16 x''[t];
dae = y[t] == y0[t];
ics = x[0] == 0, x'[0] == 0;

sol[10.] = NDSolve[ode, ics, dae, x, t, 0, 10];
xp = x' /. sol[10.];  (* iterate with next approximation to x' *)
sol["Final"] = NDSolve[ode, ics, dae, x, t, 0, 10];

Comparemos con la solución producida por el método numérico de Laplace utilizado por la respuesta de xzczd. En lo que sigue, usaremos

sol["Laplace"] = x -> FunctionInterpolation[GWR[f, t], t, $MachineEpsilon, 10]

dónde f y GWR son como en la otra respuesta.

Las soluciones son aproximadamente las mismas:

Plot[x[t] /. sol["Laplace"], x[t] /. sol["Final"], t, 0, 10]

Gráficos de Mathematica

Podemos comparar qué tan bien las soluciones siguen el ODE. La principal razón por la que el método de Laplace parece mucho peor se debe a FunctionInterpolation. Sin embargo, parece ser una mejor aproximación a valores pequeños de t. La función opODE da el residuo de una solución dada sol en el momento t de la ODE del OP, con NIntegrate en lugar de Integrate.

opODE[t_?NumericQ, sol_] := Hold[
   0.01 - 6.25 x[t] + (1.2 NIntegrate[x'[t - τ]/Sqrt[τ], τ, 0, t])/10^7 - 16 x''[t]
   ] /. sol // ReleaseHold;

GraphicsRow[
 Plot[opODE[t, sol["Laplace"]], opODE[t, sol["Final"]],
    t, ##, PlotPoints -> 20, MaxRecursion -> 2, PlotRange -> All
    ] & @@@ 0., 0.001, 0.001, 1, 1, 10
 ]

Gráficos de Mathematica

NDSolve es una mejor alternativa a FunctionInterpolation para construir una interpolación precisa. Curiosamente, el método de Laplace muestra un comportamiento errático similar cerca t == 0 como el NDSolve-método de iteración. La función GWR del paquete de inversión numérico de Laplace necesita el argumento t ser numérico, pero no lo protege con ?NumericQ; de ahí la envoltura gwr debajo. Con este método de interpolación, el método numérico de Laplace parece comparable.

gwr[t_?NumericQ] := GWR[f, t];
sol["Laplace"] = NDSolve[x[t] == gwr[t], y'[t] == 1, 
    y[$MachineEpsilon] == $MachineEpsilon, 
   x, t, $MachineEpsilon, 10];

GraphicsRow[
 Plot[opODE[t, sol["Laplace"]], opODE[t, sol["Final"]],
    t, ##, PlotPoints -> 20, MaxRecursion -> 2, PlotRange -> All
    ] & @@@ 0., 0.002, 0.002, 1, 1, 10
 ]

Gráficos de Mathematica

Presumiblemente x -> gwr (de @xzczd) produce una solución muy precisa, pero lleva mucho tiempo evaluarla. Por ejemplo,

opODE[0.01, x -> gwr] // AbsoluteTiming
opODE[9.95, x -> gwr] // AbsoluteTiming
(*
  27.7618, -1.13159*10^-13 - 7.74942*10^-15 I
  35.6744, -1.36696*10^-15 - 6.18062*10^-26 I
*)

Respuesta de Mathematica 11.0.0 for Microsoft Windows (64-bit) (July 28, 2016)

ieqn = 0.01 - 6.25*x[t] + 
1.2* Integrate[x'[t - s]/Sqrt[s], s, 0, 10]/10^7 == 16*x''[t];
ic = x[0] == 0, x'[0] == 0;
sol = DSolve[ieqn, ic, x[t], t]

$ x (t) a 0.0016 , -0.0016 cos (0.625 t) $

Plot[x[t] /. sol, t, 0, 10]

ingrese la descripción de la imagen aquí

MMA puede resolver ecuaciones integro-diferenciales ahora.


EDITADO 10.07.2018:

Parece que en MMA 11.3DSolveno puedo resolver y volver sin evaluar.Probablemente sea un error (regresión).

Te mostramos reseñas y puntuaciones

Si te sientes impulsado, puedes dejar una división acerca de qué te ha gustado de este enunciado.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *