1
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – user64494
    Apr 6 at 9:36
  • $\begingroup$ 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. $\endgroup$
    – user64494
    Apr 6 at 11:21
  • $\begingroup$ 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]. $\endgroup$
    – user64494
    Apr 6 at 16:26
3
$\begingroup$

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

Figure 1

$\endgroup$
2
  • $\begingroup$ Thank you @Alex Trounev. Is there a difference between FindRoot and NMinimize approaches? In theoretical terms, I mean. $\endgroup$
    – Mark
    Apr 6 at 8:24
  • 1
    $\begingroup$ @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. $\endgroup$ Apr 6 at 11:15
3
$\begingroup$

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).

enter image description here

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

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Apr 6 at 20:05
  • $\begingroup$ @ 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. $\endgroup$ Apr 7 at 9:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.