# How to validate the solution of a system of integro-differential equations?

I have successfully solved a system of integro-differential equations.

The system is

Eq1 = x1'[t] == (1 - 3/2*Exp[t]) x1[t] + 1/2 Exp[3*t] x2[t] +
Integrate[1/4 Exp[-t + 2*s - x1[s]^2] x1[s], {s, 0, t}];

Eq2 = x2'[t] == (-Exp[3*t]) x1[t] - 2 Exp[t] x2[t] +
Integrate[1/4 Exp[-t + 2*s - x2[s]^2] x2[s], {s, 0, t}]


Now following @J.M. comment, I convert the above system to a system of ODE's, like this,

nEq1 = x1'[t] == (1 - 3/2*Exp[t]) x1[t] + 1/2 Exp[3*t] x2[t] + Exp[-t]*x10[t]
nEq1x10 = x10'[t] == 1/4 Exp[2*t - x1[t]^2] x1[t]
nEq2 = x2'[t] == (-Exp[3*t]) x1[t] - 2 Exp[t] x2[t] + Exp[-t]*x20[t]
nEq2x20 = x20'[t] == 1/4 Exp[2*t - x2[t]^2]*x2[t]


With the initial conditions

ics = {x1 == 1, x10 == 0, x2 == 2, x20 == 0};


Now solving the converted system using NDSolve

sys = Join[{nEq1, nEq1x10, nEq2, nEq2x20}, ics];
sol = NDSolve[sys, {x1[t], x2[t]}, {t, 0, 4}];


Finally, plotting the results,

Plot[Evaluate[{x1[t], x2[t]} /. sol], {t, 0, 4}, PlotRange -> All] My question is, how we can validate the solutions?

• You want to check the validity of the solution with respect to the orginal Eq1, Eq2? How did you calculate the ODE system? Feb 5, 2017 at 16:54
• @Feyre That's right!, if possible? I converted the integral to a differential and end up with two extra odes.
– zhk
Feb 5, 2017 at 16:56
• As far as I can tell, NDSolve is quite robust for solving IVP of ODE(s), so the result will be correct as long as the deduction of ODEs is correct. BTW you can deduce x10'[t] with D[Integrate[1/4 Exp[-t + 2*s - x2[s]^2] x2[s]/Exp[-t], {s, 0, t}] , t] Feb 6, 2017 at 2:51

You can test the two sides of Eq1 and Eq2.

If the ODE is valid, the two sides of each equation are equal. You can break up the equations into parts. For the right hand sideEq1, the first two parts will be:

sol1p1 = ((1 - 3/2*Exp[t]) x1[t] /. sol) + 1/2 Exp[3*t] x2[t] /. sol


In the interim, plotting these with x1[t] becomes:

Plot[Evaluate[{D[x1[t] /. sol,
t], sol1p1}], {t, 0, 4}, PlotRange -> All] You need to add a variable replacement for the integral:

sol1p2[tt_] :=
NIntegrate[
1/4 Exp[-t + 2*s - (x1[s] /. (sol /. t -> s))^2] x1[s] /. (sol /.
t -> s), {s, 0, t}] /. t -> tt


And plotting it:

Plot[sol1p2[t], {t, 0, 4}, PlotRange -> All] Adding them together, and subtracting the original lhs of Eq1:

Plot[sol1p5[t][] + sol1p1[] - D[x1[t] /. sol, t], {t, 0, 4}] And some values:

Table[(sol1p5[t][] + sol1p1[] - D[x1[t] /. sol, t]) /.
t -> tt, {tt, 0, 4}]


{{0.}, {0.0000332613}, {-4.37529*10^-6}, {-0.0000395381}, \ {-0.0000143361}}

Which are very close to zero, as we would expect.

You can do the same for Eq2`.