I have a fairly complicated set of coupled non-linear integro-differential equations that I am trying to solve using NDSolve. The equations are:
$\dot{x}=\big|y(t)-x(t)\big|^{1/n}\left[\text{Sign}[y(t)-x(t)]-x(t)\right]$
$\big|y(t)-x(t)\big|^{1/n}\text{Sign}\left(\big|y(t)-x(t)\big|\right)=\cos(t+\pi\alpha/2)-\int\limits_0^t (t-\xi)^{\alpha-1}\ \dot{y}(\xi)\ d\xi$
with initial conditions $x(0)=0, y(0)=0$. Some of you may notice that the last integral is nothing but the definition of the fractional order integral with order $0<\alpha<1$. I implement this system of equations in Mathematica as follows
\[Alpha] = .5;
n = .5;
tmax = 10;
Intermediate[t_Real, y_Real] :=
Evaluate[Integrate[(t - \[Xi])^(\[Alpha] - 1) y, {\[Xi], 0, t},
PrincipalValue -> True]]
SolutionOfSimultaneousEquations :=
NDSolveValue[{x'[t] == (Abs[y[t] - x[t]])^(1/n)*(Sign[y[t] - x[t]] - x[t]),
(Abs[y[t] - x[t]])^(1/n) Sign[y[t] - x[t]] == Cos[t + Pi/2 \[Alpha]] - Intermediate[t, y'[t]],
x[0] == 0,y[0] == 0}, {x, y}, {t, 0, tmax}]
and then use the following bit of code to plot my results:
SolutionForx = SolutionOfSimultaneousEquations[[1]];
SolutionForxFn[t_] := SolutionForx[t];
SolutionFory = SolutionOfSimultaneousEquations[[2]];
SolutionForyFn[t_] := SolutionFory[t];
Plot[{SolutionForxFn'[t], SolutionForxFn[t]}, {t, .001, tmax},
PlotStyle -> Automatic, PlotRange -> Full]
Essentially I have defined a function called Intermediate[t_Real,y_Real]
and this function is called by NDSolve everytime the value of the integral is required. However the solution I get from this code does not make physical sense. My main concern is with the way I have implemented the integral. Notice that the integral required the value of $y(t)$ for all times between $0$ and $t$. However, currently my code only passes a single value for $y(t)$ and the integral is completely different. Is this correct?
What do I do to ensure that this integral is performed for all times $t$? I understand that only using something like an Euler forward scheme starting from $t=0$ will ensure that the integral can be performed (so that I have the value of $y$ for all prior times). Is there a better way to solve integro-differential equations?
Integrate
only integrates over a value of y'[t], thus the evaluated integral is quite different from what you want. Unfortunately, integro-differential equations are not handled out of the box byNDSolve
. $\endgroup$ – user21 May 16 '13 at 6:51data = {}; Intermediate[t_Real, y_Real] := Block[{if}, data = Join[data, {{t, y}}]; if = Interpolation[data, "InterpolationOrder" -> 1]; NIntegrate[(t - \[Xi])^(\[Alpha] - 1) if[\[Xi]], {\[Xi], 0, t}] ]
anddata = {}; NDSolve[....
However, I get anNDSolve::icfail
message. But perhaps it's a start. $\endgroup$ – user21 May 16 '13 at 16:341/0^.5
type error occurs internally triggering as a resultNDSolve::icfail
. Hope you find a way out.. $\endgroup$ – PlatoManiac May 16 '13 at 20:42