11
$\begingroup$

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: fie

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

Naive solution

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.

  • 0
    I think you might have had more answers asking this on http://stackoverflow.com (you might want to consider flagging a mod to migrate this question, though it is not off-topic here)2011-11-30

2 Answers 2

5

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 
5

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)