335 views (last 30 days)

Deal all.

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

David Goodmanson
on 28 Dec 2020

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
on 29 Dec 2020

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
on 30 Dec 2020

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
on 3 Jan 2021

Hi Emma,

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 6 Comments

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1224677

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1224677

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1225482

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1225482

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1225907

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1225907

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1225972

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1225972

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1226767

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1226767

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1229562

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/701157-please-help-me-convert-spectral-integral-equation-to-matlab-code#comment_1229562

Sign in to comment.