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

-
Yes, including the code would be a good idea. If it's too long, there's Pastebin. – J. M. Jul 27 '11 at 15:52
@J.M.: Edited, added the code with comments. – Ilya Jul 27 '11 at 16:04
I think you might have had more answers asking this on stackoverflow.com (you might want to consider flagging a mod to migrate this question, though it is not off-topic here) – Tobias Kienzler Nov 30 '11 at 8:03

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

-

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)

-