Numerical solution of an integral equation

I have problems with a solution of an integral equation in MATLAB: all conditions are double-checked, but the answer is incorrect. Let me state the equation: $$ x(s) = g(s)+\int\limits_0^1x(t)f(t-s)\,dt\quad(1) $$ where $$ f(t) = \frac{\sqrt{2}}{\pi(1+(t+1)^4)} $$ and $$ g(s) = \int\limits_{1-s}^\infty f(t)\,dt. $$ For the solution $x$ I know that it should be bounded $0\leq x(s)\leq 1$ and monotonically increasing. To solve this equation I used fie toolbox which is also supported by many publications of K. Atkinson and provides strict bounds on the error. The graph of the solution:

Note an ugly jump close to $0$ as well as non-monotonic behavior of the function. An error said to be less than $10^{-4}$.

This result didn't satisfied me and I then found $x$ with a naive approach: made a grid on $[0,1]$ with a step $10^{-3}$ and replaced $x(s)$ by $x_i = x(s_i)$, $g(s)$ by $g_i = g(s_i)$, $f(t-s)$ by $f_{ij} = f(s_j-s_i)$ and integral by a sum. As a result I obtained $(\mathbf I - \delta \mathbf f)\mathbf w = \mathbf g$. Solution of this system is on this picture. With this method an error is less than $0.01$.

I used fie toolbox for the problem which is much tough and the result was perfect, but for the current equation something goes wrong. I hope that it's due to the fact I've encoded something incorrect.

I spent already three days trying to find a mistake - but there is no so much code and I wonder if you can help me. Maybe, it will be better if I put the code here for someone to help me check it? In the case this question is off-top, I am sorry - just tell me, I'll delete it.

P.S. Mathematica is able to give an explicit form for $g$, so I can put it here if it helps.

Edited: Added Code. I use the following procedure to call fie

[soln,errest,cond] = Fie(1,0,1,1,@kPdfNew,@rhsNew)

First numbers are $\lambda =1$, $a = 0$, $b = 1$ and behavior $=1$ - smooth. For the function $g$ I use

function output = rhsNew(x)
n = length(x);
output = x;
for i=1:n
        output(i) = quad(@Poi,1-x(i),100);
end

and for the kernel I put

function output = kPdfNew(s,t)
z = t-s; 
n = length(z);
output = z;
for i = 1:n
    output(i) = Poi(z(i));
end

where Poi is a density $f(t)$ given by

function output = Poi(z)
n = length(z);
output = z;
for i=1:n
    output(i) = sqrt(2)/pi/(1+(z(i)+1)^4);
end

You may comment on the integration till $100$ rather then $\infty$ - but first, the difference $\int\limits_{100}^\infty f(t)\,dt$ is smaller then a machine precision and second, before I used the explicit form for $g$ and the result was the same.

Solution 1:

Just in the case someone will have the same problem, here is the answer. The problem was in working with vector-argument, vector-value functions in MATLAB, so the kernel should be given in a slightly different way:

function output = kPdfNew(s,t)
z = t-s; 
[nx,ny] = size(s);
output = zeros(size(s));
    for ii = 1:nx
      for jj = 1:ny
        output(ii,jj) = sqrt(2)/pi/(1+(z(ii,jj)+1)^4);
      end
    end
end

Solution 2:

Just in case performance and readability are issues, you can achieve what your answer did with the following code as well:

function output = kPdfNew(s,t)
output = sqrt(2)/pi./(1+(t-s+1).^4);
end

(The main problem with your original code was that you iterated for i=1:length(z), instead of length you could have used numel)