# Solving Fredholm Equation with composite unknown function

I would like to numerically solve a Fredholm Equation where the unknown function is composite. For example, an equation like the one described in Solving Fredholm Equation of the second kind but having composite functions as unknowns.

Consider then the Fredholm Equation: $$\phi\left(\frac{x^2}{2}-1\right) = 1 + \frac12 \int_{0}^{\pi} \text{cos}\left(x-s\right) \, \phi\left(\frac{s^2}{2}-1\right) \,ds$$ for $$x\in[0,\pi]$$.

How could one use Mathematica to find a numerical solution?

• The command of Maple 2021 restart; evalf(intsolve(phi(t) = 1 + 1/2*int(cos(sqrt(2 + 2*t) - s)*phi(s)/sqrt(2 + 2*s), s = -1 .. Pi^2/2 - 1), phi(t), method = collocation, order = 2)); produces $\phi(t)= 0.1998091304 t^{2}- 0.9895444856 t+ 0.8053645268$. The plot plot(eval(0.1998091304*t^2 - 0.9895444856*t + 0.8053645268, t = x^2/2 - 1), x = 0 .. Pi) is strikingly different from the one in the below answer. Apr 6, 2021 at 9:36
• It should be in the above: restart; evalf(intsolve(phi(t) = 1 + 1/2*int(cos(sqrt(2 + 2*t) - sqrt(2 + 2*s))*phi(s)/sqrt(2 + 2*s), s = -1 .. Pi^2/2 - 1), phi(t), method = collocation, order = 2)); produces $\phi(t)=- 0.4666943120 t^{2}+ 2.240582170 t+ 1.558351011$. The plot plot(eval(-0.4666943120*t^2 + 2.240582170*t + 1.558351011, t = x^2/2 - 1), x = 0 .. Pi) is dramatically different from the below plot. Apr 6, 2021 at 11:21
• The command of Maple 2021 restart; evalf(intsolve(phi(t) = 1 + 1/2*int(cos(sqrt(2 + 2*t) - sqrt(2 + 2*s))*phi(s)/sqrt(2 + 2*s), s = -1 .. Pi^2/2 - 1), phi(t), method = collocation, order = 3)); produces $\phi(t)= 0.1526021469 t^{3}- 1.249904734 t^{2}+ 2.400982898 t+ 3.155295017$. One may look at its plot in x by Plot[0.1526021469*t^3 - 1.249904734*t^2 + 2.400982898*t + 3.155295017 /. t -> x^2/2 - 1, {x, 0, Pi}, PlotPoints -> 30]. Apr 6, 2021 at 16:26

This equation can be solved with Bernoulli polynomials. First we substitute $$y=x^2/2-1$$, $$t=s^2/2-1$$, then define colocation points, solution and system of equations as follows

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


Numerical solution can be evaluated in two step

sol = NMinimize[Norm[eqs], A]

sol1 = FindRoot[Table[eqs[[i]] == 0, {i, Length[eqs]}],
Table[{a[i], a[i] /. sol[[2]]}, {i, Length[A]}]]


Finally we plot solution

Plot[{u[x^2/2 - 1] /. sol1, u[x^2/2 - 1] /. sol[[2]]}, {x, 0, Pi},
PlotRange -> All, Frame -> True,
PlotLegends -> {"FindRoot", "NMinimize"},
FrameLabel -> {"x", "\[Phi]"}, PlotLabel -> Row[{"n = ", nN}]]


• Thank you @Alex Trounev. Is there a difference between FindRoot and NMinimize approaches? In theoretical terms, I mean.
– Mark
Apr 6, 2021 at 8:24
• @Mark Please, pay attention, that we use sol in sol1 as initial guess. So, we can consider NMinimize as first iteration in numerical solution. The Newton's iterative method is implemented in FindRoot, while NMinimize always attempts to find a global minimum in a given space. Sometimes two solutions are not differ, but in this example there is difference between. Apr 6, 2021 at 11:15

With the transformed integral equation (thanks @Alex Trounev) collocation method(Galerkin method, doesn't need initial guess! ) with a polynomial basis up to order y^11

g=Function[y, Table[y ^i, {i, 0,  11}]
i0 = NIntegrate[ g[y] , {y, -1, Pi^2/2 - 1}  ]// Quiet;
i1 = NIntegrate[Outer[Times, g[y], g[y]], {y, -1, Pi^2/2 - 1} ]// Quiet;
i2 = 1/2 NIntegrate[Outer[Times, g[y], g[t]] Cos[ Sqrt[2 (1 + y)] - Sqrt[2 (1 + t)]]/Sqrt[2 (1 + t)], {y, -1, Pi^2/2 - 1}, {t, -1,Pi^2/2 - 1} ]// Quiet;


leads to a result

Plot[LeastSquares[i1 - i2, i0] . g[x^2/2 - 1], {x, 0, Pi},PlotRange -> {0, All}]


which shouldn't differ from Alex Trounev's answer (but slightly does).

The difference might occur because the grid in Alex Trounev's answer doesn't include the boundaries y=-1,y=Pi^2/2-1 ?

• It is nice approach (+1), but NIntegrate gives large error due to singularity at t=-1, it is why I don't use this point as a colocation point. Apr 6, 2021 at 20:05
• @ AlexTrounev Thanks for the hint! The singularity t==-1 seems to be integrable, but the condition of matrix i1-i2 isn't as bad Det[i1-i2]== O[10^25] as expected. Perhaps the function basis has to be adapted. Apr 7, 2021 at 9:19