Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
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

2 Answers 2

up vote 4 down vote accepted

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
share|improve this answer

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)

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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