Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this community
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?
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}]]
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.
$\endgroup$
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 ?
NIntegrate gives large error due to singularity at t=-1, it is why I don't use this point as a colocation point.
$\endgroup$
Apr 6, 2021 at 20:05
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.
$\endgroup$
Apr 7, 2021 at 9:19
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 plotplot(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. $\endgroup$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 plotplot(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. $\endgroup$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 inxbyPlot[0.1526021469*t^3 - 1.249904734*t^2 + 2.400982898*t + 3.155295017 /. t -> x^2/2 - 1, {x, 0, Pi}, PlotPoints -> 30]. $\endgroup$