# How to handle this constant-interval integral in ODE?

The problem is to solve a system of nonlinear integro-differential equation with a Fredholm integral. $Q_i(x), i=1,2,3, x \in [-\text{max},\text{max}]$ \begin{align*} m_i{Q_i}''(x)-V_i(\vec{Q})-d_i\int_{-\text{max}}^{\text{max}} {\frac{Q_i(x)-Q_i(x')}{(x-x')^2}dx'} = 0 \end{align*}

How can I solve the cases like $d_i=0.01,0.1$ or any other values. It appears to me Laplace transform doesn't help. I failed to turn it to a system of ODEs by differentiation, either. Or I missed something?

Update: An important check for the $d_i=0$ case is to see whether a conserved quantity $-|\vec{Q}|+\sum_i{\frac{(Q_i-w_i)^2}{\alpha_i}}-\frac{1}{2}\sum_i{m_i {Q_i}'^2}$ is constant or not.

I use the following code to easily solve the case when all $d_i=0$. The parameters are set to their typical values for my need ($\alpha_1-\alpha_{2,3}>2w_1>0=w_{2,3}$), but one can tune. $V_i=\frac{2}{\alpha_i}(Q_i-w_i)-\frac{Q_i}{|\vec{Q}|}$ and $w_i$ is constant. The form of boundary condition is not to be changed.

BTW, the integral kernal $\frac{Q_i(x)-Q_i(x')}{(x-x')^2}$ is singular, although somehow suppressed by the nominator. If not possible to solve this, adding some small cutoff $\frac{Q_i(x)-Q_i(x')}{(x-x')^2+\delta^2}$ is OK.
When all $d_i=0$, it seems that one of $Q_2,Q_3$ stays at zero. I wonder if this holds even when $d_i\neq0$.

ClearAll["Global*"];
max = 5.0; accur = 15;
w1 = 0.01; w2 = 0; w3 = 0; w = {w1, w2, w3};
α1 = 1.5; α2 = 1.0; α3 = 0.8; α = {\
α1, α2, α3};
Qstart = {w1 - α1/2, 0.0, 0.0}; Qend = {w1 + α1/2, 0.0,
0.0};
M1 = 1/α1; M2 = 1/α2; M3 = 1/α3; M = {M1, M2, M3};
dV = {2/α1 (Q1[u] - w1) - Q1[u]/Sqrt[
Q1[u]^2 + Q2[u]^2 + Q3[u]^2],
2/α2 (Q2[u] - w2) - Q2[u]/Sqrt[Q1[u]^2 + Q2[u]^2 + Q3[u]^2],
2/α3 (Q3[u] - w3) - Q3[u]/Sqrt[
Q1[u]^2 + Q2[u]^2 + Q3[u]^2]};
s = NDSolve[{M1 Q1''[u] == dV[[1]], M2 Q2''[u] == dV[[2]],
M3 Q3''[u] == dV[[3]], Q1[-max] == Qstart[[1]],
Q2[-max] == Qstart[[2]], Q3[-max] == Qstart[[3]],
Q1[max] == Qend[[1]], Q2[max] == Qend[[2]],
Q3[max] == Qend[[3]]}, {Q1, Q2, Q3}, {u, -max, max},
AccuracyGoal -> accur];
QQ = First[{Q1, Q2, Q3} /. s];
Plot[{QQ[[1]][u], QQ[[2]][u], QQ[[3]][u]}, {u, -max, max},
PlotRange -> All, PlotLegends -> Automatic]


The only solution I can think out is to make a step backward i.e. turn to finite difference method. I'll use pdetoae for the generation of difference equation.

V[i_] = (2 (Q[i][x] - w[i]))/α[i] - Q[i][x]/Sqrt@Sum[Q[i][x]^2, {i, 3}];
transrule = lst_List[i_] :> lst[[i]];
pararule =
{w -> {0.01, 0, 0}, α -> {1.5, 1, 0.8}, d -> {0.01, 0.02, 0.03},
m -> 1/α, ϵ -> 10^-6};
{Qstart, Qend} = {{w[1] - α[1]/2, 0, 0}, {w[1] + α[1]/2, 0, 0}} //.
pararule /. transrule;

Clear@kernel
kernel[Q[i_]@x_] = (Q[i][x] - Q[i][ξ])/((x - ξ)^2 + ϵ) /. pararule;

SetAttributes[int, Listable];
eq = Table[m[i] Q[i]''[x] - V[i] - d[i] int[Q[i]@x] == 0, {i, 3}] //. pararule /.
transrule;
max = 5;
points = 100;
difforder = 4;
domain = {-max, max};

{nodes, weights} = Most[NIntegrateGaussRuleData[points, MachinePrecision]];

midgrid = Rescale[nodes, {0, 1}, domain];

intrule = int[
qix_] :> -Subtract @@ domain weights.Map[Function[ξ, #] &@kernel[qix], midgrid];

grid = Flatten[{domain // First, midgrid, domain // Last}];

ptoafunc = pdetoae[{Q[1], Q[2], Q[3]}[x], grid, difforder];
del = #[[2 ;; -2]] &;
ae = del /@ ptoafunc[eq] /. intrule;
aebc = Table[
{Q[i][domain // First] == Qstart[[i]],
Q[i][domain // Last] == Qend[[i]]}, {i, 3}];


Notice that pdetoae can't handle integral so some extra coding has been made. The code for discretization of integral is modified from this answer.

Last step is to solve the nonlinear equation system with FindRoot:

sollst = With[{initial = -2},
Partition[FindRoot[{ae, aebc},
Flatten[Table[{Q[i][x], initial}, {i, 3}, {x, grid}], 1],
MaxIterations -> 200][[All, -1]], points + 2]];

ListLinePlot[{grid, #}\[Transpose] & /@ sollst, PlotRange -> All, Axes -> {False, True}]


Notice there seem to be multiple solutions for the system. With initial = -2:

With initial = 2:

• Thank you so much for this answer! It seems that you combined ingenuity and your arsenal for ODE/PDE. I myself tried combining my code with an iterative method, which turns out to be too slow to see any output. I am a novice in both differential equation and Mathematica and not sure whether I can absorb/understand your codes. I will study it with great effort, although maybe I hope to ask sometime. If possible, it would be quite helpful if you could add more explanation. Thanks again for your time and effort! – xiaohuamao Nov 18 '17 at 12:24
• It seems that the boundary conditions are satisfied well and solution is fast and stable. As for two solutions you show, first is reasonable while second is somehow trivial from my background knowledge since only one Q is varying. But the imaginary parts are not very small while the whole problem should be real. Is there a cure or some possible way to improve? And from your experience and expertise in ODEs, is this (real part) solution reliable in general? – xiaohuamao Nov 18 '17 at 12:31
• I'm not familiar with finite element method or this ODE->PDE->AE approach. I once discretized linear ODEs to solve matrix eigenvalues and hence imagine here it becomes nonlinear AE solved by FindRoot or I just misunderstood? But what's the machinery to turn it into some AE? BTW, I also encountered not small imaginary parts of eigenvalues when solving some discretized ODEs which should have real eigenvalues. But the small and large eigenvalues are usually very real if I remember correctly. A bit similar or maybe just coincidence. – xiaohuamao Nov 18 '17 at 12:38
• Sorry, am I missing something? Tuning the cutoff value, the solution doesn't change at all? – xiaohuamao Nov 18 '17 at 13:12
• @huotuichang Oh it's actually not a "ODE->PDE->AE" approach, it's just a differential equation -> difference equation i.e. finite difference method (FDM) approach (though I've named my function as pdetoae). You just need to remember pdetoae is something that'll replace derivative terms with corresponding difference formulas. If there's really a real solution, maybe we need a better initial guess for FindRoot, but I just tried to use the solutions when $d_i=0$ as initial guesses (yeah if you adjust options of "Shooting" in NDSolve, you'll be able to find 2 different solutions)… – xzczd Nov 18 '17 at 13:12