Please help me convert spectral integral equation to matlab code.

Deal all.

I need you help to convert this spectral integral equation to matlab code

[Walter Roberson]

What kind of integral is that? I think I recently saw that notation for something having to do with Laplace transforms ?

[David Goodmanson]

Hello Emma,

This appears to be a one-dimensional integral in nu, in which case the e^(-gamma1*|x2-y2|) factor can be taken out of the integral as a constant. That doesn't seem quite right. What is the domain of integration, L?

[Walter Roberson]

L is the portion of the complex plain in which the absolute value of the complex component is less than some constant ?

[Walter Roberson]

I would expect a circular region for that formula, not a rectangular "strip", since x and y are complex valued -- the Re{v} and the Im{v} and the "Complex v - plane" on the diagram tell us that.

[David Goodmanson]

Hi Emma,

Could you provide more information on the integral? For example, is the integral an indefinite integral, or an integral along the entire line from -inf-i*b to inf-i*b? Also, there appear to be either branch points or poles in the nu plane, the effects of which are not evident in the integrand, and it's not so clear what x and y are. On the face of it, if x and y are constant then the integral is easy, so I confess that I don't have a real picture yet of what is being integrated.

[David Goodmanson]

Hello Emma,

Everything through

d2 G_hat / d x2^2 -gamma^2 G_hat = -exp(-i*nu*y1)*delta(x2-y2)

looks good.

Suppose the integral listed in the question is called A. You have written k^2(x2). Is it true that k^2 is a function of x2? If so, then A, with its factor of e^(-gamma1*abs(x2-y2)), is not the correct solution. Instead you will have to solve a second order differential equation, which might not be easy.

If k^2 is constant, then the solution is usually expressed in terms of bessel functions, with no explicit integration required. But if you want to do this by integration, then it looks like the A integral is basically going to work. But I'm not convinced that the four branch points you show are valid. It would be best to proceed and construct A from the ground up to establish the branch points, but first I wanted to check and see if k^2 is indeed a constant.

Answer 1

[David Goodmanson]

Hi Emma,

Seeing the actual problem helped.

It sounds like this is supposed to be a region-1-only problem with region 2 out of the picture. Assuming that’s the case, you can find k^2, and k itself. Of course k is complex, as is gamma1 = sqrt(nu.^2-k^2).

It’s possible to use the ‘integral’ function to do the integration, but that function has some problems because the branch point involving 1/gamma1 is close to the axis of integration and causes sharp peaks. It works to use the simple trapz integration, and there is more control over how it’s done.

Once you have k, here is a sample calculation

In the problem, they use (x1,x2) for the source point coordinates and (y1,y2) for the observation point coordinates. This is really confusing. It should have been (x1,y1) for source and (x2,y2) for obs. Well, we are stuck with the system.

nu = -50:.00001:50;
gamma1 = sqrt(nu.^2-k^2);
x1 = 1;
x2 = 1;
y1 = 3;
y2 = 2;
f = (1/(4*pi))*exp(i*nu*(x1-y1)).*exp(-gamma1*abs(x2-y2))./gamma1;
figure(1)
plot(nu,real(f),nu,imag(f))
grid on
 
T = trapz(nu,f)
B = (i/4)*besselh(0,1,k*sqrt((x1-y1)^2+(x2-y2)^2))

The plot allows you to see the extent of the integrand f, and to pick upper and lower limits on nu to make sure that f is tiny there. You can see large peaks in f, and to successfully integrate across those takes small spacing in nu, such as the .00001 here. You can reduce the spacing until the answer T does not change significantly.

If you compare the resulting integral T to the bessel expression B, they agree to about 12 significant figures in this case, which is not bad.

[David Goodmanson]

Hi Emma,

can't answer that without knowing more about the original question, because the bessel expression is not there.

Possibility A: You are in region 1 only, neglecting the fact that there is a region 2. In this case the bessel expression will be something very simple,

(i/4)*besselh(0,1,k*sqrt((x1-y1)^2+(x2-y2)^2))

or something close to it, where k is in region 1.

Possibility Z: You have true Sommerfeld scattering, in which case the bessel expression will be complicated and will involve k from both region 1 and region 2.

If it's A, then since y2>0 and x2>0 (terrible notation) the problem seems to be telling you that both source and obs are in region 1. In that case all you have to do is plug gamma1 (with k from region 1) into the expression for f and use trapz. The trapz calculation agrees with the bessel calculation, as you can check. That's it. The integration only gives you one (x1,x2,y1,y2) case at a time, so you would have to do a lot of cases to build up enough points to do a plot. But I don't think the problem is asking for that. It's just asking for comparison at a few points.

[David Goodmanson]

Hi Emma,

All I can say is, you just need to try something. In my code I was in region 1, I calculated k_1^2 and gamma1, did the integral and it agreed with the bessel expression. That pretends that region 2 does not exist, but at least you get an anwer.

You expression above is doing the same thing for region 2 and pretending that region 1 does not exist. You can find k_2^2 in region 2, which is provided in the problem, find gamma2, plug that into the code and see that it agrees with the bessel expression. Again, you get an answer.

However, it looks like the instructor wants a solution assuming both regions 1 and 2 are present. I do not how to to that and it may well be difficult.

[Osman Senturk]

Hi Emma,

It looks you removed all your comments. Actually i I just wanted to some help.