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]
Observación:
-
eq
esRationalize
d porque el primer argumento deGWR
no debe contener números aproximados. -
$MachineEpsilon
se elige como el límite inferior deFunctionInterpolation
porqueGWR
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
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]
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
]
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
]
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]
MMA puede resolver ecuaciones integro-diferenciales ahora.
EDITADO 10.07.2018:
Parece que en MMA 11.3
DSolve
no 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.