# Numerical solution of Fredholm Equation

I would like to plot the solution of the Fredholm Equation $$f\left(x\right)+\frac{1}{\pi} \int_{-1}^{1} \frac{1}{1+\left(x-t\right)^2}f\left(t\right) dt=1, \ \ (|x|\leq 1)$$ I tried to use Mathematica to find a numerical solution:

PHI = DSolveValue[\[Phi][x] == 1 - 1/Pi*Integrate[\[Phi][t]/(1 + (x - t)^2),
{t, -1, 1}], \[Phi], x]
Plot[PHI, {x, -1, 1}]


But I obtained some errors and I can't plot the numerical solution. How can I fix this problem?

• In your Latex line there is one equation and in the Mathematica code there is second one. What actually do you try to solve? Apr 4 '21 at 23:03
• Your code works for me for the DSolveValue part, but you won't get a plot without assigning a value for lambda. If you mean what you say in your code, you need to change your Latex to match. Apr 5 '21 at 0:40
• Sorry everyone, I was wrong to copy the function in Mathematica. I corrected the post.
– Mark
Apr 5 '21 at 6:34
• Are you looking for analytical solution only? This equation can be solved by numerical method like colocation or iteration. Apr 5 '21 at 9:50
• Just to compare. The command of Maple 2021 intsolve(f(x) + int(f(t)/(1 + (x - t)^2), t = -1 .. 1)/Pi = 1, f(x), method = collocation, order = 3) performs $$f \! \left(x \right) = 0.09775266381 x^{2}+ 0.6581514078 .$$ Analytical output is too long to be cited here. Apr 5 '21 at 12:20

Perhaps NestList gives an iterated solution (Picard iteration)

sol = NestList[1 - 1/Pi*Integrate[#/(1 + (x - t)^2), {t, -1, 1},Assumptions -> -1 <= x <= 1] &, 0 , 10];

Show[ Plot[sol  , {x, -1, 1}, PlotRange -> {0, 1}],
Plot[sol[[-1]], {x, -1, 1}, PlotStyle -> {Thickness[.005], Black}]]


This gives the symbolic Neumann-Series .

• As we know, for some Fredholm equations the iteration method is not converge at n->Infinity. It has minimal error at n about 6. What do you think about this equation? Apr 5 '21 at 12:48
• @ AlexTrounev How did you check the error? If I compare the difference of two succesive iterations the numerical error seems to degree monoton. Apr 5 '21 at 12:58
• The error is the difference between the true solution and an iterative one, not between two iterations. Apr 5 '21 at 14:17
• @AlexTrounev But the true solution isn't known so far? Apr 5 '21 at 14:34
• Yes it is. The true solution is known in theory only. Apr 5 '21 at 15:51

You have a wonderfully behaving smooth almost flat kernel. Since you want a numerical solution, then why not simply solve it with a quadrature?

# To matrix equation with quadrature

Here's a sample code in WolframAlpha notebook:

n=100;
K=Table[1/(1 + (2/n*(i - j))^2), {i,n}, {j, n}];
A=N[IdentityMatrix[n]+2/n*1/Pi*K];
f=LinearSolve[A,ConstantArray[1,n]];
ListPlot[N[Table[{2/n*(i -1/2)-1,f[[i]]}, {i,n}]]]


Here's how it works. Convert the integral equation into a matrix equation replacing the integral with the dumbest rectangle quadrature: $$\int K(x,t) f(t)dt=\sum_{i=1}^NK(x,t_i)f(t_i)\frac 2 N$$ where $$t_i=2/N*(i-1/2)-1$$.

Now, after plugging the same knots into $$x$$ variable you get a linear system: $$f_{t_i}+\frac 2 N \sum_jK_{t_it_j}f_{t_j}= \sum_j\left(\delta_{ij}+\frac 2 N K_{t_it_j}\right)f_{t_j}=1$$ $$\left(I+\frac 2 N K\right)f=\mathbf 1$$ which can be solved trivially with high precision as shown above.

You can go fancy by employing Gaussian quadrature. Here's an example:

n = 27
<< NumericalDifferentialEquationAnalysis;
knot = 0.5
gw = Join[gwL,gwC,gwR];
n = Length[gw]

K=Table[1/(1 + ((gw[[i,1]] -gw[[j,1]]))^2)*gw[[j,2]], {i,n}, {j, n}];
A=N[IdentityMatrix[n]+1/Pi*K];
f=LinearSolve[A,ConstantArray[1,n]];
ListPlot[N[Table[{gw[[i,1]],f[[i]]}, {i,n}]]]

f[[Floor[n/2]+1]]

Out:     0.657412


# Spectral decomposition

There is a different solution that gives an insight into the integral kernel. Apply the spectral decomposition to the kernel: $$-\frac 1 \pi \frac 1 {1+(x-t)^2}=\sum_{i=1}^\infty\lambda_i \psi_i(x)\psi_i(t)$$ Once you get the eigen values $$\lambda_i$$ and vectors $$\psi_i(x)$$, the solution to your problem is: $$f(x)\approx\sum_{i=1}^K \frac 1 {1-\lambda_i}\left(\int_{-1}^1\psi_i(t)dt\right)\psi(x)$$ , where $$K$$ is the number of eigenvectors used for approximation.

Usually, we order eigenvalues descending by absolute value. Because the kernel is smooth and almost flat only the first few eigenvalues are required because they decay so quickly:

The first few eigen vectors are interesting to observe. The first one gets very close to the solution:

... and clearly only odd ones are symmetric and should have nonzero values of $$\int_{-1}^1\psi_i(t)dt$$ integral as shown next. Also, these integrals quickly decay too, helping to cut $$K$$, the number of required eigen vectors for high precision approximation:

We can see now that if the right hand side wasn't simply a constant, 1 in your case, then a different - potentially larger - number $$K$$ of eigenvectors may have been required for a good approximation.

Here's the full code:

n=999;
x=N[Table[2*(i -1/2)/n-1, {i,n}]];
d=N[Table[1/(1 +(2/n* (i -j))^2)/Pi,{i,n},{j,n}]];
K = -d;
A=N[2/n*K];
{vals,v} = Eigensystem[A,11];

v1 = v[[1]]*Sum[v[[1,i]],{i,n}]/(1-vals[[1]]);
v3 = v[[3]]*Sum[v[[3,i]],{i,n}]/(1-vals[[3]]);
v5 = v[[5]]*Sum[v[[5,i]],{i,n}]/(1-vals[[5]]);
v7 = v[[7]]*Sum[v[[7,i]],{i,n}]/(1-vals[[7]]);
v9 = v[[9]]*Sum[v[[9,i]],{i,n}]/(1-vals[[9]]);
v11 = v[[11]]*Sum[v[[11,i]],{i,n}]/(1-vals[[11]]);
vapp = v1+v3+v5+v7+v9+v11;
ListPlot[Table[{x[[i]],vapp[[i]]},{i,n}]]

ListLogPlot[-vals]
ListPlot[Table[{x[[i]],v[[1,i]]}, {i,n}]]
ListPlot[Table[{x[[i]],v[[2,i]]}, {i,n}]]
ListPlot[Table[{x[[i]],v[[3,i]]}, {i,n}]]
ListPlot[Table[{x[[i]],v[[4,i]]}, {i,n}]]


One more colocation method based on Bernoulli (or Euler) polynomials with using NMinimize and FindRoot

 nN = 12; xcol = Table[-1 + 2 j/nN, {j, 0, nN}];
v[x_] := Table[BernoulliB[n, x], {n, 0, nN}]; A = Array[a, {nN + 1}];
u[x_] = A . v[x];
eqs = Table[
u[xcol[[i]]] - 1 +
1/Pi A .
NIntegrate[v[t]/(1 + (xcol[[i]] - t)^2), {t, -1, 1},
AccuracyGoal -> 8], {i, Length[xcol]}];


First method

sol = NMinimize[Norm[eqs], A]


Second method

sol1 = FindRoot[Table[eqs[[i]] == 0, {i, Length[eqs]}],
Table[{a[i], 1/10}, {i, Length[A]}]]


Visualization

Plot[{u[x] /. sol1, u[x] /. sol[[2]]}, {x, -1, 1}]


• For x=0 both methods yield 0.65741153600362792,but using method user @Ulrich Neumann give us 2/3.Which is correct one? Apr 5 '21 at 10:20
• @MariuszIwaniuk Actually 2/3=0.66666..., and colocation method gives solution with some error of $h^2$. In my code h=1/12, so we have error of 1/144. Apr 5 '21 at 10:57
• I tried with nN = 36 and give me: 0.6574115... ? Apr 5 '21 at 11:30
• @MariuszIwaniuk Ah, I see, that your method also gives 0.657412 for n=10. May be true value is not 2/3? Apr 5 '21 at 12:07
• Looks like at x=0 value is not 2/3.Try: N[Expand[DSolveValue[\[Phi][x] == 1 - 1/Pi* Integrate[\[Phi][ t]*(Series[1/(1 + (x - t)^2), {x, 0, 4}] // Normal) // Simplify // Expand, {t, -1, 1}], \[Phi][x], x]], 20] Apr 5 '21 at 15:50

Long ago I borrowed the code from site, but I don't know exactly where.

ClearAll["*"]; Remove["*"];

n = 5;
a1 = -1;
b1 = 1;
B[x_, i_] := Binomial[n, i]*((x - a1)^i*(b1 - x)^(n - i))/(b1 - a1)^n;
k[t_, x_] := -(1/(1 + (-t + x)^2));
a[x_] := 1;
f[x_] := 1;
\[Lambda] = 1/Pi;

nIntegrate[most__] := NIntegrate[most, Method -> {Automatic, SymbolicProcessing -> 0}];

K2[j_?NumericQ] := nIntegrate[B[x, j] f[x], {x, a1, b1}] // Quiet;
mid[x_?NumericQ, i_?NumericQ] := nIntegrate[\[Lambda] k[t, x] B[t, i], {t, a1, b1}] // Quiet;

K3[j_?NumericQ, i_?NumericQ] :=
nIntegrate[(a[x] B[x, i] - mid[x, i]) B[x, j], {x, a1, b1}] // Quiet;

sol = Table[Total@Table[M[i] K3[j, i], {i, 0, n}] == K2[j], {j, 0, n}];

func = Total@Table[M[i] B[x, i], {i, 0, n}] /. NSolve[sol, Table[M[i], {i, 0, n}]][[1]]

Plot[func, {x, -1, 1}, PlotRange -> {Automatic, {0, 1}}]
`

• Some hours ago there was a tricky numerical solution at the end of your answer , why did you delete it? Apr 5 '21 at 13:05
• @Ulrich Neumann I'm not sure it was correct. Apr 5 '21 at 13:59
• What a pity I was hoping it would be something like the "Method of lines for DAE systems"! What was wrong? Apr 5 '21 at 14:03