# Solving a Volterra integral equation numerically

I would like to solve for $P(t)$, in Mathematica, a Volterra integral equation of the 2nd kind. It is:

$$P(t) = R_0(t) + \int_0^t P(t') R_0(t-t')dt'$$

I know the function $R_0$ and would like to find a general way to solve for $P(t)$ given any $R_0$.

I am aware that it is a convolution and can be solved in Laplace domain but I do not want to do it that way since inverting the transform brings about a whole other set of problems for me.

• Volterra equations (1st and 2nd kind) are a common solution to many electrochemical problems. I have solutions for these in my book "Simulating Electrochemical Reactions with Mathematica" but the short answer is to adopt a finite difference scheme using what in the echem literature is known as Huber's method. The solution can be obtained very rapidly using standard approaches in Mma. – Mike Honeychurch Apr 24 '12 at 4:42
• @MikeHoneychurch I looked at your book an Amazon, very interesting! Maybe you can give a brief outline as an answer (wouldn't want to plagiarize your book, of course). Although this question lacks specificity, the general techniques (or references) might be of interest... – Jens Apr 24 '12 at 5:36
• Here is Mathematica Journal article about solving integral equations. Perhaps that is useful. – user21 Apr 24 '12 at 7:31
• @Jens I went and had a look at the relevant parts of the book prior to posting. Because of the way the thing is laid out it didn't seem like something that could readily be chopped up and displayed here to readily answer the Ops question. If the OP is familiar with discretizing an integral equation then he should be able to work through the problem. Next step is to code it. I'll see if I can work through some generic examples. – Mike Honeychurch Apr 24 '12 at 7:33
• The purpose of what I wrote back then was a step by step discussion of how to implement Hubers method. I never put it together into a function or anything like that. Would be happy to email the OP the relevant bits. – Mike Honeychurch Apr 24 '12 at 7:44

Mathematica is an incredible tool for checking conjectures and making sketches. I'm going to demonstrate it below.

Let's start with checking that in case when $R_0$ (I replaced it with $R$) is a polynomial, the solution of this Volterra equation reduces to linear ODE. Lets take some derivatives of the equation:

ClearAll[P, R, s, t];
eqn = P[t] == R[t] + Integrate[P[s] R[t - s], {s, 0, t}]
Table[D[eqn, {t, n}], {n, 0, 2}] // TableForm


We see that this process is somehow similar to integration by parts. If $R^{(n)}$ is zero here, we get an ODE. Let's check it out.

R = 1 + 2 #^2 &;
deg = Exponent[R[t], t];
Table[D[eqn, {t, n}], {n, 0, deg + 1}] // TableForm


The required initial conditions are obtained using intermediate derivatives:

iconds = Table[D[eqn, {t, n}] /. t -> 0, {n, 0, deg}];
TableForm[ndeqs = {D[eqn, {t, deg + 1}]}~Join~iconds]


Now we can solve this either numerically or symbolically:

dsol = DSolve[ndeqs, P, t][[1, 1]]
ndsol = NDSolve[ndeqs, P, {t, 0, 1}][[1, 1]]
GraphicsRow[Plot[P[t] /. #, {t, 0, 1}, PlotRange -> All] & /@ {dsol, ndsol}]


Verifying that the symbolic solution is exact:

eqn /. dsol // Simplify
(* ==> True *)


Now, if kernel is not polynomial we can just interpolate it taking Chebyshev points of second kind for higher accuracy:

RR = Exp[-#] + #*Sin[#] - #*Cos[#^2] &;
deg = 9;
nodes = Table[(1 + Cos[(j \[Pi])/deg])/2, {j, 0, deg}] // N;
R = Evaluate[InterpolatingPolynomial[Transpose@{nodes, RR /@ nodes}, #]] &;
Plot[{RR[t] - R[t]}, {t, 0, 1}, PlotStyle -> Thickness[Large]]


As we see the error is rather small. Here I took kernel from Daniel Lichtbau's answer. Now we can reuse the above code and obtain the following plot of approximate solution:

This perfectly agrees with Daniel's results. Thank you.

Something similar showed up on StackExchange a while back (last August).

https://stackoverflow.com/questions/6974929/how-can-i-reference-a-specific-point-of-my-function-inside-ndsolve

Here is a variant from one response that I cobbled together to illustrate your case. I chose the kernel function (your R_0) to be e^(-x) + x sin(x) - x cos(x^2). I attempt to solve the integral equation over the range from 0 to 1, iterating to get improvements. We start by setting the initial value of P(x) equal to kernel(x). We'll do an iteration by evaluation the integral with this initial P(x), at 100 points in (0,1). We use those to make an interpolating function, and that becomes the new approximation to P(x). I did four iterations of this process.

Clear[ifunc]
kernel[x_] := Exp[-x] + x*Sin[x] - x*Cos[x^2]
func[x_, 0] := kernel[x]
ifunc[0][x_] := kernel[x]

func[x_?NumericQ, n_Integer] :=
kernel[x] +
NIntegrate[kernel[x - y]*ifunc[n - 1][y], {y, 0, x},
MinRecursion -> 2, AccuracyGoal -> 3]

ifunc[j_Integer /; j >= 1] := ifunc[j] = Module[{vals},
vals = Table[{x, func[x, j]}, {x, 0, 1, .01}];
Interpolation[vals]]

t = AbsoluteTime[];
Table[Print[{j, AbsoluteTime[] - t}];
t = AbsoluteTime[];
ifunc[j], {j, 0, 4}];

{0,0.000479}

{1,0.000019}

{2,0.995606}

{3,26.250936}

{4,26.335029}


A plot indicates that the last two are indistinguishable. That does not by any means imply that what I did is correct, but simply that it converged. To something.

Plot[Evaluate[Table[ifunc[j][x], {j, 0, 4}]], {x, 0, 1},
PlotStyle -> {Blue, Red, Green, Purple, Black}]


There are many sources that give practical advice on how to program such an integral equation. So I went a different route and asked myself how to get a solution by using the most literal application of the defining equation.

By that I mean no explicit discretization, keeping the integral. That can be done using the assumption that the fixed-point theorem holds for the given equation. That theorem means that if I replace $P(t')$ under the integral by the entire expression for $P(t)$, with $t\to t'$, and keep doing that with all the $P(...)$ that appear, the result will converge.

Be warned that this is best suited for the short-time limit.

For two iterations, that would look like this:

Clear[r]; i = 0; FixedPoint[(i++; (r[K[i + 1]] +
Integrate[# r[K[i] - K[i + 1]], {K[i], 0, K[i + 1]}])) &,
r[K[1]], 2]


$r(K[3])+\int_0^{K[3]} r(K[2]-K[3]) \left(r(K[2])+\int_0^{K[2]} r(K[1]) r(K[1]-K[2]) \, dK[1]\right) \, dK[2]$

Here, I ran FixedPoint for only 2 steps because it's a symbolic calculation (I'll make it numerical soon, that's why I went with FixedPoint). The K[i] are symbolic integration variables.

Below is a randomly chosen kernel:

r[t_] := t^2


It's not at all realistic, but the integrals are easy to do:

i = 0; FixedPoint[(i++; (r[K[i + 1]] +
Integrate[# r[K[i] - K[i + 1]], {K[i], 0, K[i + 1]}])) &,
r[K[1]], 4]


$\frac{K[5]^{14}}{2724321600}+\frac{K[5]^{11}}{2494800}+ \frac{K[5]^8}{5040}+\frac{K[5]^5}{30}+K[5]^2$

This is the desired approximation for $P(t)=P(K[5])$. Here, K[5] is the time variable (the outermost integration limit).

Next, I'll make a numerical function out of this, so we can in principle deal with more complicated r[t]:

Clear[step];
step[t_?NumericQ, function_] :=
r[t] + NIntegrate[function[\[Tau]] r[t - \[Tau]], {\[Tau], 0, t}]


This defines one iteration step for the FixedPoint search below:

solution[time_] :=
FixedPoint[Function[{t}, step[t, #]] &, r,
SameTest -> (Abs[#1[time] - #2[time]] < 10^(-3) &)][time]


To limit the time of execution, I set a SameTest that isn't too stringent about when it considers convergence to be achieved.

The result is as follows:

l = Table[solution[t], {t, 0, 2, .1}]


{0., 0.0100003, 0.0400107, 0.090081, 0.160341, 0.251042, 0.362595, 0.495614, 0.650956, 0.829768, 1.03353, 1.26411, 1.5238, 1.81539, 2.14222, 2.50824, 2.91812, 3.37726, 3.89198, 4.46953, 5.11828}

This is only a proof of principle, not intended to be of practical use... I just posted it because it's a concise formulation.