MATLAB Answers

Please help me convert spectral integral equation to matlab code.

335 views (last 30 days)
Emma Rice
Emma Rice on 23 Dec 2020
Commented: Osman Senturk on 3 Jan 2021
Deal all.
I need you help to convert this spectral integral equation to matlab code

  6 Comments

Show 3 older comments
Walter Roberson
Walter Roberson on 24 Dec 2020
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
David Goodmanson on 24 Dec 2020
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
David Goodmanson on 27 Dec 2020
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.

Sign in to comment.

Answers (1)

David Goodmanson
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.

  3 Comments

David Goodmanson
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
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
Osman Senturk on 3 Jan 2021
Hi Emma,
It looks you removed all your comments. Actually i I just wanted to some help.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by