Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: integral equation with singularities and double integral
Date: Mon, 31 Oct 2011 13:31:28 +0000 (UTC)
Organization: Saab Underwater Systems AB
Lines: 35
Message-ID: <j8m7vg$1m$1@newscl01ah.mathworks.com>
References: <j8c3fp$8kv$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-02-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1320067888 54 172.30.248.47 (31 Oct 2011 13:31:28 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Mon, 31 Oct 2011 13:31:28 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 920395
Xref: news.mathworks.com comp.soft-sys.matlab:747704

"Sulaiman Rabbaa" wrote in message <j8c3fp$8kv$1@newscl01ah.mathworks.com>...
> Hello
> I need help to solve the following integral equation:
> double integral (with respect to x and y) of 137.03.*y.^2./((0.238.*exp(0.067.*y.^2)+1).*(w-5.26.*x.*y-2.63).*(w-5.26.*x.*y+2.63))=1+8478./(10828-w.^2-1.13.*j.*w)
> xmin=-1, xmax.=1, ymin=0, ymax=inf 
> I want to find w which is complex number.
> I tried the following code:
> f=@(w) dblquad(@(x,y) 137.03.*y.^2./((0.238.*exp(0.067.*y.^2)+1).*(w-5.26.*y.*x-2.63).*(w-5.26.*y.*x+2.63)),-1,1,0,100)-1-8478./(10828-w.^2-j.*w.*1.13)
> V=fsolve(f,160)
> There are some singularities. How can I solve the problem

Hi,
The Newton-Raphson solution below gives w = (78.1093 - 0.2275i). Another solution is w = (-78.1093 - 0.2275i).

f = @(x,y,w) 137.03.*y.^2./(0.238.*exp(0.067.*y.^2)+1) ...
    ./(w-5.26.*x.*y-2.63)./(w-5.26.*x.*y+2.63);

df = @(x,y,w) 137.03.*y.^2./(0.238.*exp(0.067.*y.^2)+1) ...
    ./(w-5.26.*x.*y-2.63).^2./(w-5.26.*x.*y+2.63).^2 ...
    .*2.*(5.26.*x.*y-w);

g = @(w) quad2d(@(x,y)f(x,y,w),-1,1,0,100);
dg = @(w) quad2d(@(x,y)df(x,y,w),-1,1,0,100);

h = @(w) g(w) - 1 - 8478./(10828-w.^2-1.13i.*w);
dh = @(w) dg(w) - 8478.*(2*w+1.13i)./(10828-w.^2-1.13i.*w).^2;

w = 80-1i;
dw = 1;
while abs(dw) > 1e-9 && abs(w) < 1e6
    dw = -h(w)/dh(w);
    w = w + dw;
end

/Jonas