I am trying to find a numerical solution for an equation of the form:

$$ f(t) = \int_{t_{min}}^{t} \mathrm{d}t' {\exp[2 (t^\prime - t)] E(t^\prime) f(t^\prime)} + \int_{t}^{0}{\mathrm{d}t' \exp(t^\prime - t) E(t^\prime) f(t^\prime)} $$

where $E$ is the complete elliptic integral of the first kind and $t_{min}$ can be in principle any negative integer. I have tried NIntegrate and some iteration but didn't get me too far... can anyone help me with a better idea?

Thank you very much for your help!!!

  • $\begingroup$ you might want to look at mathematica.stackexchange.com/a/11609/1089 mathematica.stackexchange.com/questions/9544/… $\endgroup$
    – chris
    Commented Dec 7, 2012 at 6:19
  • $\begingroup$ @chris, I have already seen the links, but it seems that the standard methods to solve Fredholm equations somehow don't work in this case :(. can you be a little more explicit with your second suggestion? Thank you very much! $\endgroup$
    – mia
    Commented Dec 7, 2012 at 6:32
  • $\begingroup$ @chris, I have tried your suggestion but the problem is that the function is very badly oscillating around $t=0$ and I am not really able to keep it under control. If you have time and if I am not asking too much, would you like try it out yourself and tell me how it looks like, in particular for small values of t? Thank you! $\endgroup$
    – mia
    Commented Dec 7, 2012 at 7:33
  • $\begingroup$ Also possibly relevant is this stack overflow post $\endgroup$ Commented Dec 7, 2012 at 15:23
  • $\begingroup$ The integral equation can be turned into the 2nd order DE, $f''(t)+3f'(t)+(2+E(t))f(t)=0$, (modulo mistakes) but the initial/boundary conditions are a little tricky to work out. $\endgroup$
    – Simon
    Commented Dec 7, 2012 at 22:27

2 Answers 2


This might help you get started. I think this is a variation on a method called Frobenius' method. The idea is to use the fact that the equation is linear to expand the solution over a set of basis functions (here B-Splines) and find the corresponding coefficients. It should provide you with an approximate solution (which you can improve upon while adding more Splines in your basis), provided the sought solution is smooth enough and the Kernel is regular enough.

Let us define some sampling for the sought solution

np = 5; tmin = -5; Δt = -tmin;
kfun[n_, d_] := Join[ConstantArray[0, d], Range[0, 1, 1/(n - d)],  ConstantArray[1, d]];
knots = tmin + Δt*kfun[np + 1, 3];

Let us format the unknown $F(t)$ represented by a sum of B-Splines

Format[a[i_]] = Subscript[a, i];
F[t_] = Sum[BSplineBasis[{3, knots}, i, t] a[i], {i, 0, np}]

and define the corresponding basis

basis[t_] = Table[BSplineBasis[{3, knots}, i, t], {i, 0, np}]

Let's show the basis just because we can:

Plot[basis[t] // Evaluate, {t, -5, 0}]

Mathematica graphics

Now your equation is (in terms of F[t]=f[t] EllipticK[t])

eqn = F[t]/EllipticK[t] == Integrate[Exp[2 (t - s)] F[s], {s, tmin, t}] +
   Integrate[Exp[(t - s)] f[s], {s, t, 0}]

It involves the convolution kernel

Kern[s_, t_] = 
 Piecewise[{{Exp[2 (s - t)], s < t}, {Exp[(s - t)], s > t}}]

Now let us define the matrix of dot products of our basis over the kernel K

M = 
   basis[t][[i + 1]] basis[s][[j + 1]] Kern[s, t], {t, tmin, 0}, {s, 
    tmin, 0}],
  {i, 0,np}, {j, 0,np}]

and the matrix of dot products of our basis (as it is not ortho-normal)

Q = ParallelTable[
   basis[t][[i + 1]] basis[t][[j + 1]]/EllipticK[t] , {t, tmin, 0}],
  {i, 0,np}, {j, 0,np}]

Now let us look at the eigen-space of that equation

{eig, vec} = 
  Inverse[Q].M - IdentityMatrix[Length[Q]] // Eigensystem;

Our approximate solution is corresponds to the eigenvector with the smallest eigenvalue:

var = Table[a[i], {i, 0, np}]; ra = Thread[var -> Last[vec]];

Plot[F[t] /EllipticK[t] /. ra // Evaluate, {t, -5, 0}]

Mathematica graphics

Note that this solution is defined up to a (possibly negative) multiplicative constant.

  • 1
    $\begingroup$ Thank you so much Chris! I am trying to continue along your idea but I am quite slow -- mathematica is not really my forte, I have resorted to it after having failed in fortran :(. It is very frustrating, as I am stuck in this equation since a while now... Again, thank you so very much for your help, I hope I will manage to get it done... $\endgroup$
    – mia
    Commented Dec 7, 2012 at 11:08
  • $\begingroup$ Sorry for annoying you again, but it seems to be a problem with the null space, which is always empty regardless of the kernel (I have tried several). Do you know what is happening? $\endgroup$
    – mia
    Commented Dec 7, 2012 at 14:53
  • $\begingroup$ I have also corrected some bugs in your code, still doesn't do, but I am very tired so maybe I just don't see things. Please let me know if you notice anything and again, thank you for taking so much time... $\endgroup$
    – mia
    Commented Dec 7, 2012 at 18:09
  • $\begingroup$ Note that I was a bit optimistic to assume that the approximate solution (the sum of BSplines) was exactly in the NullSpace of Q-M. Now I just look to the closest null solution. $\endgroup$
    – chris
    Commented Dec 7, 2012 at 19:40
  • $\begingroup$ Chris, I know this is not the right place to say thank you, but apparently mia is not able to post anymore (no idea why), so I am just going to say that your code works just beautifully! I had spotted a couple of things that needed to be corrected, but you did it before I have noticed the possibility of editing somebody else's answer :). Thank you! ps. this is another mia, I just opened another account -- for unknown reasons the old mia is not able to write comments anymore :) $\endgroup$
    – mia
    Commented Dec 7, 2012 at 21:40

The commenters got it almost solved, missing only that $t_{\text{min}}$ is an unknown instead of a parameter. The question should be: for what values of $t_{\text{min}}$ does the integral equation have a nontrivial solution?

As shown in the comments, the function satisfies

$$ [2 + E(t)] f(t) + 3 f'(t) + f''(t) = 0\\ 2 f(0) + f'(0) = 0\\ f(t_{\text{min}}) + f'(t_{\text{min}}) = 0 $$

All of these equations are linear, so the solution can be multiplied by an arbitrary scalar and still remains a solution. This last observation allows us to set $f(0)=1$ as an additional boundary condition. Together with the first two conditions above, we thus get the solution

sol = NDSolve[{(2 + EllipticE[t]) f[t] + 3 f'[t] + f''[t] == 0, 
               2 f[0] + f'[0] == 0, f[0] == 1}, f, {t, -10, 0}];
F = f /. First[sol];
Plot[F[t], {t, -10, 0}]

enter image description here

The allowed values for $t_{\text{min}}$ are those where this solution satisfies the third condition $f(t_{\text{min}}) + f'(t_{\text{min}}) = 0$:

Plot[F[tmin] + F'[tmin], {tmin, -10, 0}]

enter image description here

which we can find numerically:

FindRoot[F[tmin] + F'[tmin] == 0, {tmin, -0.7}]
(*    {tmin -> -0.657447}    *)

FindRoot[F[tmin] + F'[tmin] == 0, {tmin, -3}]
(*    {tmin -> -2.92468}    *)

FindRoot[F[tmin] + F'[tmin] == 0, {tmin, -5}]
(*    {tmin -> -4.95128}    *)

These are the "eigenvalues" of the given integral equation and represent the set of discrete values of $t_{\text{min}}$ for which the integral equation has a nontrivial solution. There are infinitely many of them, stretching towards $-\infty$.

Numerical Stability & Automatic Solving

The function $f(t)$ oscillates with exponentially increasing amplitude. To regularize the problem, we can introduce the function

$$ g(t) = e^{\frac32t}f(t) $$

which satisfies

$$ [E(t)-1/4] g(t) + g''(t)= 0\\ g(0) + 2g'(0) = 0\\ g(t_{\text{min}}) - 2g'(t_{\text{min}}) = 0 $$

Solving for $g(t)$ gives a regularly oscillating function:

sol = NDSolve[{(EllipticE[t] - 1/4) g[t] + g''[t] == 0,
               g[0] + 2 g'[0] == 0, g[0] == 1}, g, {t, -20, 0}];
G = g /. First[sol];
Plot[G[t], {t, -20, 0}]

enter image description here

As @J.M. points out, WhenEvent can be used to extract the allowed values for $t_{\text{min}}$ directly,

sol = NDSolve[{(EllipticE[t] - 1/4) g[t] + g''[t] == 0, 
               g[0] + 2 g'[0] == 0, g[0] == 1, 
               WhenEvent[g[t] - 2 g'[t] == 0, Print[t]]},
              g, {t, -20, 0}];

(*    -0.657447
      -19.7327    *)

Use a Sow/Reap combination instead of Print in order to use these values later in a calculation.

  • 1
    $\begingroup$ Of course, one can use WhenEvent[] to get a pile of admissible values for $t_{\mathrm{min}}$. $\endgroup$ Commented May 23, 2020 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.