Numerical solution of Fredholm Equation

Question

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?

Answer 1

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 .

Answer 2

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.

Gauss-Legendre quadrature

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

n = 27
<< NumericalDifferentialEquationAnalysis`;
knot = 0.5
gwL = GaussianQuadratureWeights[n/3, -1, -knot];
gwC = GaussianQuadratureWeights[n/3, -knot, knot];
gwR = GaussianQuadratureWeights[n/3, knot,1];
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}]]

Answer 3

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}]

Answer 4

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}}]