Group Abstract

Message Boards

WOLFRAM COMMUNITY

4.7K Views

7 Replies

1 Total Like

View groups...

Follow this post

Share this post:

GROUPS:

Mathematics

How to calculate the numerical integration that contains singular point?

Henan Fang

Posted 6 years ago

G is a numerical integration,the codes are as following. I found that the curve of the output is not smooth and seems like noise. This is because the integrand Abs[A2]^2 has singular point. How to solve this problem? Many thanks! Clear["`"] vh = 1; d1 = 1.3410^(-9); d2 = 1.3410^(-9); mu = 5.5; HBAR = 1.0545726610^(-34); ME = 9.109389710^(-31); ELEC = 1.6021773310^(-19); Kh = 2.9510^(10); kc = Sqrt[2MEELEC/HBAR^2]; k := kcSqrt[mu] kh := Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2] khg := Sqrt[ k^2 - (2KhSin[afa/2]Sin[afa/2] - k Sin[x] Cos[y])^2 - (2KhSin[afa/2]Cos[afa/2] + k Sin[x] Sin[y])^2] khpl := Sqrt[ k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + kc^2 vh] khplpl := Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + 2kc^2 vh] khgplpl := Sqrt[k^2 - (2KhSin[afa/2]Sin[afa/2] - k Sin[x] Cos[y])^2 - (2KhSin[afa/2]Cos[afa/2] + k Sin[x] Sin[y])^2 + 2kc^2 vh] A2 := Exp[ Ikhpld1]/(Exp[I(kh + khgplpl - khg)d1] + Exp[Ikhplpld1]) G := Re[NIntegrate[ k Sin[x] Abs[A2]^2, {x, 0, Pi/2}, {y, -Pi/6, Pi/6}]]; tmr := {afa, G}; Export["D://1.34.txt", Table[tmr, {afa, 0, Pi/6, 0.005}], "Table"];

G is a numerical integration,the codes are as following. I found that the curve of the output is not smooth and seems like noise. This is because the integrand Abs[A2]^2 has singular point. How to solve this problem? Many thanks!

Clear["`*"]
vh = 1;
d1 = 1.34*10^(-9);
d2 = 1.34*10^(-9);
mu = 5.5;
HBAR = 1.05457266*10^(-34);
ME = 9.1093897*10^(-31);
ELEC = 1.60217733*10^(-19);
Kh = 2.95*10^(10);
kc = Sqrt[2*ME*ELEC/HBAR^2];
k := kc*Sqrt[mu]
kh := Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2]
khg := Sqrt[
  k^2 - (2*Kh*Sin[afa/2]*Sin[afa/2] - 
      k Sin[x] Cos[y])^2 - (2*Kh*Sin[afa/2]*Cos[afa/2] + 
      k Sin[x] Sin[y])^2]
khpl := Sqrt[
  k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + kc^2 vh]
khplpl := 
 Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + 
   2*kc^2 vh]
khgplpl := 
 Sqrt[k^2 - (2*Kh*Sin[afa/2]*Sin[afa/2] - 
      k Sin[x] Cos[y])^2 - (2*Kh*Sin[afa/2]*Cos[afa/2] + 
      k Sin[x] Sin[y])^2 + 2*kc^2 vh]
A2 := Exp[
   I*khpl*d1]/(Exp[I*(kh + khgplpl - khg)*d1] + Exp[I*khplpl*d1])
G := Re[NIntegrate[
    k Sin[x] Abs[A2]^2, {x, 0, Pi/2}, {y, -Pi/6, Pi/6}]];
tmr := {afa, G};
Export["D://1.34.txt", Table[tmr, {afa, 0, Pi/6, 0.005}], "Table"];

POSTED BY: Henan Fang

7 Replies

Sort By:

Hans Dolhaine

Hans Dolhaine, retired

Posted 6 years ago

Your code without numbers gives (without integration) Clear["`"] kc = Sqrt[2MEELEC/HBAR^2]; k := kcSqrt[mu] kh := Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2] khg := Sqrt[ k^2 - (2KhSin[afa/2]Sin[afa/2] - k Sin[x] Cos[y])^2 - (2KhSin[afa/2]Cos[afa/2] + k Sin[x] Sin[y])^2] khpl := Sqrt[ k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + kc^2 vh] khplpl := Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + 2kc^2 vh] khgplpl := Sqrt[k^2 - (2KhSin[afa/2]Sin[afa/2] - k Sin[x] Cos[y])^2 - (2KhSin[afa/2]Cos[afa/2] + k Sin[x] Sin[y])^2 + 2kc^2 vh] A2 := Exp[ Ikhpld1]/(Exp[I(kh + khgplpl - khg)d1] + Exp[Ikhplpld1]) int = k Sin[x] Abs[A2]^2; res = DeleteDuplicates[Cases[int, Sqrt[x__ ], Infinity]] Perhaps it is an idea to check the single terms for their order of magnitude. And perhaps you may want to check res /. Sqrt[x_] :> Sqrt[FullSimplify[Expand[x]]]

Your code without numbers gives (without integration)

Clear["`*"]
kc = Sqrt[2*ME*ELEC/HBAR^2];
k := kc*Sqrt[mu]
kh := Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2]
khg := Sqrt[
  k^2 - (2*Kh*Sin[afa/2]*Sin[afa/2] - 
      k Sin[x] Cos[y])^2 - (2*Kh*Sin[afa/2]*Cos[afa/2] + 
      k Sin[x] Sin[y])^2]
khpl := Sqrt[
  k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + kc^2 vh]
khplpl := 
 Sqrt[k^2 - (Kh - k Sin[x] Cos[y])^2 - k^2 Sin[x]^2 Sin[y]^2 + 
   2*kc^2 vh]
khgplpl := 
 Sqrt[k^2 - (2*Kh*Sin[afa/2]*Sin[afa/2] - 
      k Sin[x] Cos[y])^2 - (2*Kh*Sin[afa/2]*Cos[afa/2] + 
      k Sin[x] Sin[y])^2 + 2*kc^2 vh]
A2 := Exp[
   I*khpl*d1]/(Exp[I*(kh + khgplpl - khg)*d1] + Exp[I*khplpl*d1])
int = k Sin[x] Abs[A2]^2;
res = DeleteDuplicates[Cases[int, Sqrt[x__ ], Infinity]]

Perhaps it is an idea to check the single terms for their order of magnitude.

And perhaps you may want to check

res /. Sqrt[x_] :> Sqrt[FullSimplify[Expand[x]]]

POSTED BY: Hans Dolhaine

Henan Fang

Posted 6 years ago

Thanks very much for your reply. I have tried to improve the precision of the calculation. Unfortunately, it does not work.

POSTED BY: Henan Fang

Hans Dolhaine

Hans Dolhaine, retired

Posted 6 years ago

Look at your parameters, you are using values between 10^(-34) and 10^10 HBAR = 1.0545726610^(-34); Kh = 2.9510^(10) i.e. 44 orders of magnitude. Could it be that you are having trouble with rounding errors due to insufficient precision in your calculations?

POSTED BY: Hans Dolhaine

Gianluca Gorni

Gianluca Gorni, University of Udine

Posted 6 years ago

You may try a change of integration variables that makes the singularity simpler. Your function depends only on `Sin[x] Cos[y]` and `Sin[x] Sin[y]`, if I am not mistaken. You may take these as new variables. I have not tried it myself.

POSTED BY: Gianluca Gorni

Henan Fang

Posted 6 years ago

Thanks very much for your suggestion!

POSTED BY: Henan Fang

Bill Nelson

Posted 6 years ago

I am not convinced that Abs is the issue. Look at afa=4/20; (all your other code down to but not including the NIntegrate and last two lines) Plot3D[Re[k Sin[x] Abs[A2]^2], {x, 0, Pi/2}, {y, -Pi/6, Pi/6}] and you can substitute other values up to about 1/2 for afa and see similar results. It looks like you only have one denominator in your code and afa appears repeatedly in that denominator. I expect that denominator going to zero is responsible for this behavior. With some experimentation, can you verify this is the source of the problem? Knowing that, is there anything you can do about that denominator going to zero?

POSTED BY: Bill Nelson

Henan Fang

Posted 6 years ago

Thanks very much for your reply. Yes, the meaning of the singular points in my original discussion just means that the denominator going to zero. However, even the denominator may go to zero for some points, the integration indeed converges, and the curve with varying afa should not seem like noise. I want to know what can I do to solve this problem. Many thanks!

POSTED BY: Henan Fang

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback