# Numerically solve an integro-differential equation

So the problem is to numerically solve this Integro-Differential Equation, $$-v \frac{\:\mathrm{d}\kern0.1mm x(\zeta)}{\:\mathrm{d}\kern0.1mm \zeta} + \psi(\zeta) x(\zeta) + \int_{0}^{\zeta} f(\zeta, \eta)x(\eta)\:\mathrm{d}\kern0.1mm \eta + g(\zeta)x(1) = 0 \\$$ where $v$ and $x(1)$ are known constants, $x(0)$ is an unknown constant I need to find, and the functions, $\psi(\zeta), f(\zeta, \eta), g(\zeta)$ are chosen by me. From what I have found in the site so far, these are what I have tried.

1. One of the best solutions was the direct, easy method to solve such equations in v11. However, this only worked in the simplest case when I let $v = A_0 =\psi(\zeta) = f(\zeta, \eta) = g(\zeta) =1$ and the problem was very simple. In any more complicated case, it fails. Hence, I need a numerical solver.
2. Trying to do the same thing above with NDSolveValue instead gave a bunch of errors (can show later if needed)
3. The solution here is, unfortunately, using LaplaceTransform and fails
4. The solution here is good but perhaps the functions are too complex to take the LaplaceTransform or InverseLaplaceTransform, even numerically, giving a another bunch of errors.

FYI, I am currently trying it with $$v = 1, x(1) = 1, \psi(\zeta) = \zeta, f(\zeta, \eta) = \zeta + \eta, g(\zeta) = \zeta^2$$

• @zhk Hmm, sorry, I am not sure how to take the derivative of the integrand when deriving with respect to $\zeta$ (I changed the f now; was wrong earlier). Jun 17 '17 at 3:45
• @xzczd Not really sure how they're using those hash tags in defining the equations and if it's possible to write the integro differential in that way. Jun 17 '17 at 3:45
• You don't need to change the form of equation, notice the posts linked by me and @zhk essentially points to the same method: differentiate your equation to make it a "pure" differential equation. "I am not sure how to take the derivative of the integrand" Just differentiate it directly, with D. Jun 17 '17 at 6:08
• The Volterra integral question is not about an integro-differential equation, and the validate question is about an integro-differential equation that can be easily converted into an ODE. I think this question should be reopened because it is different from both of the referenced question/answers. In particular, I don't think it's possible to convert the OP integro-differential equation into an ODE. Plus, I have a suggested solution that doesn't apply to the other questions. Jun 17 '17 at 17:18
• @CarlWoll Can you email me your answer if it's not reopened? I am also not sure how they're using the intermediate derivatives to find the required ICs in the Volterra question. I think the MATLAB idsolver could solve the problem but unfortunately there are no examples and Mma turns out to be easier in general for such things. Jun 17 '17 at 19:18

This integro-differential equation can be solved with the method mentioned in this answer i.e. differentiate the equation to make it a pure ODE.

First, interprete the equations to Mathematica code. (BTW, if you had given the Mathematica code form of the equation in your question, your question would have attracted more attention. )

v = 1; ψ[ζ_] = ζ; f[ζ_, η_] = ζ + η; g[ζ_] = ζ^2;

bc[0] = x[1] == 1;

eq = -v D[x[ζ], ζ] + ψ[ζ] x[ζ] +
Integrate[f[ζ, η] x[η], {η, 0, ζ}] + g[ζ] x[1] == 0 /. Rule @@ bc[0]
(* ζ^2 + Integrate[(ζ + η)*x[η], {η, 0, ζ}] + ζ*x[ζ] - Derivative[1][x][ζ] == 0 *)

Then, differentiate the equation, twice.

neweq = D[eq, ζ, ζ]
(* 2 + 3*x[ζ] + 2*Derivative[1][x][ζ] + 2*ζ*Derivative[1][x][ζ] +
ζ*Derivative[2][x][ζ] - Derivative[3][x][ζ] == 0 *)

neweq is a 3rd order ODE, while currently we only have 1 b.c. i.e. bc[0], so we need to deduce another two b.c.s from eq. This can be easily achieved by setting ζ to 0 in eq and D[eq, ζ]:

bc[1] = eq /. ζ -> 0
(* -Derivative[1][x][0] == 0 *)
bc@2 = D[eq, ζ] /. ζ -> 0
(* x[0] - Derivative[2][x][0] == 0 *)

Finally, solve the equation and find $x(0)$:

sol = NDSolveValue[{neweq, bc /@ Range[0, 2]}, x, {ζ, 0, 1}]
sol[0]
(* 0.232727 *)

# Update 1: Work-around for $f(\zeta ,\eta )=e^{(\zeta +1) \eta }$

As pointed out by Carl Woll, when the form of $f$ becomes more complicated, it may be impossible to differentiate the equation to a pure ODE. Still, there exists a work-around at least for $f(\zeta ,\eta )=e^{(\zeta +1) \eta }$, that is, approximating $f$ with its series expansion.

f[ζ_, η_] = Exp[(ζ + 1) η];

nmax = 10;
approx[ζ_, η_] = Normal@Series[f[ζ, η], {ζ, 0, nmax}];
(* Error Check: *)
Plot3D[f[ζ, x] - approx[ζ, x] // Evaluate, {x, 0, 1}, {ζ, 0, 1}, PlotRange -> All]

eq = -v D[x[ζ], ζ] + ψ[ζ] x[ζ] +
Integrate[approx[ζ, η] x[η], {η, 0, ζ}] + g[ζ] x[1] == 0 /. Rule @@ bc[0];

neweq = D[eq, {ζ, nmax + 1}];

bclst = Table[D[eq, {ζ, n}] /. ζ -> 0, {n, 0, nmax}];
sol = NDSolveValue[{neweq, bc@0, bclst}, x, {ζ, 0, 1}, WorkingPrecision -> 16];
sol[0]
(* 0.1498546695665442 *)

# Update 2: A simpler and probably more general work-around

It turns out to be unnecessary to approximate $f$ with Taylor series, we can differentiate the original equation enough times and then simply take away the remaining Integrate[……]:

order = 6;

eq = -v D[x[ζ], ζ] + ψ[ζ] x[ζ] +
Integrate[f[ζ, η] x[η], {η, 0, ζ}] + g[ζ] x[1] ==
0 /. Rule @@ bc[0];

rule = HoldPattern@Integrate[__] :> 0;

neweq = D[eq, {ζ, order + 1}] /. rule;

bclst = Table[D[eq, {ζ, n}] /. ζ -> 0 /. rule, {n, 0, order}];

solnew = NDSolveValue[{neweq, bc@0, bclst}, x, {ζ, 0, 1},
WorkingPrecision -> 16]; // AbsoluteTiming
(* {1.017286, Null} *)
solnew[0]
(* 0.1498546688616941 *)

But why this works? Isn't it just a coincidence?

No, it's not.

A qualitative explanation is, every time the equation is differentiated, the Integrate[……] term will play a less important role in the new equation.

A quantitative explanation is, if we approximate $f$ with a n-th order piecewise polynomial (e.g. the polynomial here), then after differentiate the equation for n+1 times, the Integrate[…] term will exactly be 0.

Though I haven't tested it, I believe this approach is more general than the Taylor expansion approach, because when the form of $f$ becomes even more complicated (e.g. piecewise) Taylor expansion may not be suitable.

Finally, an error check:

help[ζ_?NumericQ, sol_] :=
NIntegrate[E^((1 + ζ) η) sol[η], {η, 0, ζ},
Method -> {Automatic, "SymbolicProcessing" -> 0}]

test[sol_] :=
Subtract @@ eq /. HoldPattern@Integrate[__] -> help[ζ, sol] /. x -> sol
Plot[#, {ζ, 0, 1}] & /@ test /@ {solnew, cur[8]} // GraphicsRow

As one can see, this approach is more accurate than the iterative method given by Carl.

• Spamming to get the link to move this discussion to chat since it turns out I need to tailor it to MATLAB and I didn't understand a bunch of stuff in the answer, though I am sure it works, as you showed. Jun 20 '17 at 20:52
• @ZackFair …You can ask in the comment before the chat link appears, maybe we can finish it before moving to the chat. Jun 21 '17 at 5:22
• OK, how do you get the new BCs? You mention "This can be easily achieved by setting ζ to 0 in eq and D[eq, ζ]" but why is that so? Is it some constraint on the system or an assumption? Jun 23 '17 at 1:23
• @ZackFair It's just a mathematical deduction. Notice when $ζ=0$ the integration becomes zero. If you still have difficulty in understanding this step, try doing it by hand. Jun 23 '17 at 1:38
• Oh yes then it's just a derivative term and its variables in the equation, which would be a BC at $\zeta=0$. Thanks. That makes sense. Jun 24 '17 at 22:05

If you can convert the integro-differential equation into an IVP or BVP ODE, then that would be the best approach. If you can't do so, then you can try the following iterative approach. The basic idea is to replace the unknown $x(\eta )$ in the integrand with an ansatz (guess), and then solve the resulting differential equation. Repeat, with the solution to the ODE as the next guess to be used.

Here are the suggested parameters:

ν = 1;
ψ[ζ_] := ζ
f[ζ_, η_] := ζ + η
g[ζ_] := ζ^2

Given a guess cur[n-1], the next guess would be:

reset[] := (
Clear[cur];
cur[n_] := cur[n] = NDSolveValue[
{
-ν x'[ζ] + ψ[ζ] x[ζ] + Integrate[f[ζ,η] cur[n-1][η], {η, 0, ζ}] + g[ζ] == 0,
x[1] == 1
},
x,
{ζ, 0, 1}
]
)
reset[]

Starting with an initial guess:

cur[0] = #&;

The following plots show the differences in the iterates, and the final answer:

Plot[Evaluate@Table[cur[n][t]-cur[n-1][t], {n, 1, 8}], {t, 0, 1}, PlotLegends->Range[8]]
Plot[cur[8][t], {t,0,1}]

With the choice of $f(\zeta ,\eta )=\zeta +\eta$, the integro-differential equation can be turned into an ODE, as shown by @xzczd, so let's compare:

sol = NDSolveValue[
{
2 + 3 x[ζ] + 2 x'[ζ] + 2 ζ x'[ζ] + ζ x''[ζ] - x'''[ζ] == 0,
x[1] == 1,
x'[0] == 0,
x[0] - x''[0] == 0
},
x,
{ζ, 0, 1}
];

Here is the solution of the ODE:

Plot[sol[t], {t, 0, 1}]

And a plot of the difference:

Plot[sol[t]-cur[8][t],{t,0,1}]

As an example of an integro-differential equation that can't be turned into an ODE (or at least, I don't think it can), consider:

f[ζ_, η_] := Exp[(ζ+1) η]
reset[] (* to clear out old solution *)
cur[0] = #&;
Plot[Evaluate@Table[cur[n][t]-cur[n-1][t], {n, 1, 7}], {t, 0, 1}, PlotLegends->Range[7]]
Plot[cur[7][t], {t,0,1}]

The iterative solution to this integro-differential equation converges quite a bit more slowly.

• Thanks for this answer. I think it's more comprehensive since it deals with both transformable and non-transformable cases (xzczd's code didn't work for your 2nd example). I was just wondering if, in your second code block, you forgot the x[1] with the g[ζ] in the line, -ν x'[ζ] + ψ[ζ] x[ζ] + Integrate[f[ζ,η] cur[n-1][η], {η, 0, ζ}] + g[ζ] == 0, Jun 18 '17 at 16:03
• Also, in your case, cur[8][0] or cur[7][0] gives the value of $x(0)$, right? Jun 18 '17 at 16:22
• I think including x[1] in the ODE is an issue, because Mathematica is expecting all args of a function to be the same. If you want to vary x[1], you can add the variable x1 and use x[1] == x1. Jun 18 '17 at 16:41
• Yes, cut[n][0] is the nth estimate of x[0] Jun 18 '17 at 16:45
• I am not having $x(1)$ as a variable but in my original equation, the last term is $g(\zeta)x(1)$, not just $g(\zeta)$, which is what the code above seems to assume. Or am I misunderstanding? Jun 18 '17 at 16:57