# Solving IDE system using FDM based on NDSolveFiniteDifferenceDerivative

I wanted to solve the following system of integro-differential equations (IDEs) using finite difference method (FDM).

$$y_1''(t)+t^2y_1(t)-y_2''(t)+\int\limits_0^t[(t-x)y_1(x)+y_2(x)]\mathrm{d}x=(2+t^2)\mathrm{e}^t-t-\cos t+\sin t$$

$$4t^3y_1'(t)+6t^2y_1(t)+y_2'''(t)+\int\limits_0^t[y_1(x)+(t+x)y_2(x)]\mathrm{d}x=\sin t-(1+2t)\cos t+\mathrm{e}^t(1+6t^2+4t^3)+t-1$$

with initial conditions $$y_1(0)=y_1'(0)=1$$, $$y_2(0)=y_2''(0)=0$$ and $$y_2'(0)=1$$.

The exact solutions to this system are $$y_1(t)=\mathrm{e}^t$$ and $$y_2(t)=\sin t$$.

I have used the following code developed by xzczd in another post :

 max = 2;

SetAttributes[{int1, int2}, Listable];
eq = {Derivative[y1][t] + t^2*y1[t] - Derivative[y2][t] +
int1[t] - ((2 + t^2)*E^t - t - Cos[t] + Sin[t]),
4*t^3*Derivative[y1][t] + 6*t^2*y1[t] + Derivative[y2][t] +
int2[t] - (Sin[t] - (1 + 2*t)*Cos[t] +
E^t*(1 + 6*t^2 + 4*t^3) + t - 1)} == 0;
kernel11[t_, x_] = (t - x)*y1[x] + y2[x];
kernel21[t_, x_] = y1[x] + (t + x)*y2[x];
bc = {y1 == 1, Derivative[y1] == 1, y2 == 0,
Derivative[y2] == 1, Derivative[y2] == 0};

points = 25;
difforder = 5;
domain = {0, max};
{nodes, weights} =
Most[NIntegrateGaussRuleData[points, MachinePrecision]];
midgrid = Rescale[nodes, {0, 1}, domain];
intrule1 =
int1[t_] :> (-Subtract @@ domain)*
weights . (kernel11[t, #1] & ) /@ midgrid;
intrule2 =
int2[t_] :> (-Subtract @@ domain)*
weights . (kernel21[t, #1] & ) /@ midgrid;

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

ptoafunc = pdetoae[{y1[t], y2[t]}, grid, difforder];
fullae = ptoafunc[eq] /. Flatten[{intrule1, intrule2}];
aebc = ptoafunc[bc];
{blst, mat} =
CoefficientArrays[Flatten[{fullae, aebc}],
Flatten[{y1 /@ grid, y2 /@ grid}]];
sollst = LeastSquares[N[mat], -blst];

sol1 = Interpolation[Transpose[{grid, sollst[[1 ;; Length[grid]]]}]];
sol2 = Interpolation[
Transpose[{grid, sollst[[Length[grid] + 1 ;; 2*Length[grid]]]}]];

Plot[{E^y1, Re[sol1[y1]]}, {y1, 0, max}, PlotRange -> {Full, {-5, 5}}]
Plot[{Sin[y2], Re[sol2[y2]]}, {y2, 0, max}, PlotRange -> All]


where pdetoae[] can be found here. After ploting the exact functions and the FDM solutions I found they are not matching at all.

$$\mathrm{e}^t$$ $$\sin t$$ The orange coloured plots are the solutions from FDM. The plots are for $$t\in[0,2]$$.

I believe I could not write the code correctly for kernel integration, thus, seeking any kind of help especially from xzczd since the user has developed this excellent subroutine.

Modified code:

Tried with the following code as suggested by xzczd

 int[expr_, {t_, L_, R_, step_}] :=
step*Total[Table[expr, {t, L + step, R, step}]]

step = 1/10;
bL = 0; bR = 2;
grid = Table[i, {i, bL, bR, step}];

eq = {Derivative[y1][t] + t^2*y1[t] -
Derivative[y2][t] - ((2 + t^2)*E^t - t - Cos[t] + Sin[t]),
4*t^3*Derivative[y1][t] + 6*t^2*y1[t] +
Derivative[y2][
t] - (Sin[t] - (1 + 2*t)*Cos[t] + E^t*(1 + 6*t^2 + 4*t^3) + t -
1)} == 0;

kernel11 = int[(t - x)*y1[x] + y2[x], {x, 0, t}];
kernel21 = int[y1[x] + (t + x)*y2[x], {x, 0, t}];

bc = {y1 == 1, Derivative[y1] == 1, y2 == 0,
Derivative[y2] == 1, Derivative[y2] == 0};

kernelSet11 =
Transpose[{Table[
kernel11, {t, bL, bR, step}] /. {x, bL, a_} :> {x, bL, a,
step}}];
kernelSet21 =
Transpose[{Table[
kernel21, {t, bL, bR, step}] /. {x, bL, a_} :> {x, bL, a,
step}}];

difforder = 4;
ptoafunc = pdetoae[{y1[t], y2[t]}, grid, difforder];
fullae = ptoafunc[eq] +
Transpose[ArrayFlatten[{{kernelSet11, kernelSet21}}]];
aebc = ptoafunc[bc];

{blst, mat} =
CoefficientArrays[Flatten[{fullae, aebc}],
Flatten[{y1 /@ grid, y2 /@ grid}]];
sollst = LeastSquares[N[mat], -blst];
sol1 = Interpolation[Transpose[{grid, sollst[[1 ;; Length[grid]]]}]];
sol2 = Interpolation[
Transpose[{grid, sollst[[Length[grid] + 1 ;; 2*Length[grid]]]}]];

Plot[{E^y1, Re[sol1[y1]]}, {y1, 0, bR}, PlotRange -> All]
Plot[{Sin[y2], Re[sol2[y2]]}, {y2, 0, bR}, PlotRange -> All]


But it failed also. Must be doing something wrong.

• 1. The integral cannot be handled in this way because its upper limit is $t$ rather than a constant, if you need to use FDM for this problem, one possible reference is this post. 2. If FDM is not necessary, just differentiate the equation system twice. 3. The analytic solution for this problem is not Exp and Sin: test = Block[{int1, int2}, int1[t_] := Integrate[kernel11[t, x], {x, 0, t}]; int2[t_] := Integrate[kernel21[t, x], {x, 0, t}]; eq]; test /. y1 -> Exp /. y2 -> Sin // Simplify – xzczd Aug 22 '19 at 7:34
• @xzczd Sorry there were few mistakes in the second equation, which I have corrected now. I think the exact solutions will now satisfy the right hand sides. Thanks for the reference of variable upper limit. However, I have another query. What if the kernel limits are constants (say $a$ and $b$) but the independent variable $t$ varies from $0$ to $\infty$? – Soumyajit Roy Aug 22 '19 at 13:33
• That's not a problem. First solve the problem in $[a,b]$ using the method above. Then for $[b, \infty)$, the integral becomes a constant, and it's easy to handle. – xzczd Aug 22 '19 at 14:01
• @xzczd I tried with right Riemann sum but could not produce correct results. New code has been added to the question. Sorry for seeking suggestion at every step. – Soumyajit Roy Aug 23 '19 at 5:45
• fullae = ptoafunc[eq] + Transpose[ArrayFlatten[{{kernelSet11, kernelSet21}}]];  This is apparently incorrect, don't forget you've defined eq as something like {…, …} == 0. Also, those Transpose, ArrayFlatten, etc. are redundant. – xzczd Aug 23 '19 at 6:32

You're almost there, I think the main trouble you're facing is that you're not yet familiar enough with list manipulation of Mathematica. Check carefully about how I modify the code:

int[expr_, {t_, L_, R_, step_}] := step Total[Table[expr, {t, L + step, R, step}]]

step = 1/100;
bL = 0; bR = 2;
grid = Table[i, {i, bL, bR, step}];

eq = {y1''[t] + t^2 y1[t] - y2''[t] - ((2 + t^2) E^t - t - Cos[t] + Sin[t]),
4 t^3 y1'[t] + 6 t^2 y1[t] +
y2'''[t] - (Sin[t] - (1 + 2 t) Cos[t] + E^t (1 + 6 t^2 + 4 t^3) + t - 1)};

kernel11 = int[(t - x) y1[x] + y2[x], {x, 0, t}];
kernel21 = int[y1[x] + (t + x) y2[x], {x, 0, t}];

bc = {y1 == 1, y1' == 1, y2 == 0, y2' == 1, y2'' == 0};

kernelSet11 = Table[kernel11, {t, bL, bR, step}] /. {x, bL, a_} :> {x, bL, a, step};
kernelSet21 = Table[kernel21, {t, bL, bR, step}] /. {x, bL, a_} :> {x, bL, a, step};
difforder = 4;
ptoafunc = pdetoae[{y1[t], y2[t]}, grid, difforder];
fullae = ptoafunc[eq] + {kernelSet11, kernelSet21};
aebc = ptoafunc[bc];

{blst, mat} =
CoefficientArrays[Flatten[{fullae, aebc}], Flatten[{y1 /@ grid, y2 /@ grid}]];
sollst = LeastSquares[N[mat], -blst];

{sol1, sol2} = ListInterpolation[#, grid] & /@ Partition[sollst, Length@grid];

Plot[{E^y1, sol1[y1]}, {y1, 0, bR}, PlotRange -> All,
PlotStyle -> {Automatic, {Red, Dashed}}]
Plot[{Sin[y2], sol2[y2]}, {y2, 0, bR}, PlotRange -> All,
PlotStyle -> {Automatic, {Red, Dashed}}] • Yeah...that's great. I really face difficulties with list manipulation. Thanks again for making the necessary correction in the code. I am also not very familiar with FDM. Started reading the material suggested by you. – Soumyajit Roy Aug 23 '19 at 14:34