Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this community
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.
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).
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.
