8
$\begingroup$

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?

Improve this question
$\endgroup$
5
  • $\begingroup$ Note, that in your equation you integrate over y'[\Xi], however as you correctly mention the 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 by NDSolve. $\endgroup$ – user21 May 16 '13 at 6:51
  • $\begingroup$ Yes, I see that this is the problem. What I really need is a numerical scheme that starts from the initial conditions, uses some chosen value of step size $h$, and incrementally calculates the values of the variables at each subsequent step. That way, to find the value of $x(t+nh)$ and $y(t+nh)$ all we need is $x(t+(n-1)h)$ and $y(t+(n-1)h)$. I think we call this an Euler-Forward scheme? Are there any good implementations of such a scheme that I can just adapt for my specific case? $\endgroup$ – G. H. Hardly May 16 '13 at 15:21
  • 1
    $\begingroup$ Here is a wild thought: It may be possible to construct and interpolation function from 0 until t and integrate over that: data = {}; 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}] ] and data = {}; NDSolve[.... However, I get an NDSolve::icfail message. But perhaps it's a start. $\endgroup$ – user21 May 16 '13 at 16:34
  • $\begingroup$ @ruebenko Actually in this way too the 1/0^.5 type error occurs internally triggering as a result NDSolve::icfail. Hope you find a way out.. $\endgroup$ – PlatoManiac May 16 '13 at 20:42
  • $\begingroup$ Hi All, Thanks very much for your suggestions. I think the best possible route for me to take is to write my own Euler-Forward scheme where I write $f_{k+1}=f_k+h\dot{f_k}$ and so on. Yes, I will first have to figure out how to reduce my system of equations to that sort of explicit form by differentiating them, etc. I think I can manage my own Euler-Forward scheme, but if not, I will be back here for help! $\endgroup$ – G. H. Hardly May 16 '13 at 21:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.