Message Boards Message Boards

0
|
11466 Views
|
6 Replies
|
1 Total Likes
View groups...
Share
Share this post:

Numerical integration

salangnew
Posted 11 years ago
I tried to make numerical integration of function for continuous (7.nb) and discrete times (9.nb). I have the following questions concerning discrete  time:
1.  with  logarithmic scale values of time more than 4?? are not shown, but I need to get less than 1???.
2. with linear scale,  on the contrary there appears ???????? 2us, although in the original array the maximum value is 1us
3.  while calculating coefficient  (at the bottom of file) the program says about dixision by 0. But it is not possible, as in Mathcad I get the result of about 10^-12
4. function z is not calculated, but it exists analytically

http://zalil.ru/34911397/5efb62f5.52dc96a8/9.nb
9.nb
 c = 3*10^8;H = 7.2*10^5;\[Tau] = 47*10^-6;Fd = 7200;Tp = 50*10^-6;f = 3.5*10^8;
 v = 7500;Q = 18000;L = 100;s2 = 0.003;s1 = 2; \[CapitalTheta]0 = 0.018;m = 256;
 \[Lambda] = 0.022;n = 3;lm = 25;d = 1.2;\[Lambda] = 0.022;n = 3;P0 = 25;\[Eta] = 0.8;
 KIP = 0.4;kf = Sqrt[0.5];\[CapitalTheta]s = 1;
 t = Sort@{-10^-8, -20^-9, -10^-9, 0, 10^-11, 10^-10, 10^-9, 10^-8,     20*10^-7, 2*10^-6, 10^-6};
 p = Sum[NIntegrate[ Exp[(\[Pi] Fd^2 Tp^2 (m - Abs[k])^2 x^2)/H^2 - (n^2 \[Pi]^2 d^2 y^2)/
     (\[CapitalTheta]0^2 \[Lambda]^2 \H^2*141^2) - (100 lm^2 x^2 s1^2)/(H^2 L^2 \[Pi]^2 s2^2) -
      (100 lm^2 \y^2*s1^2)/(H^2 L^2 \[Pi]^2 s2^2) - (5.55 (y^2 + x^2))/(H^2 \[CapitalTheta]0^2) -
      \[Pi] (f^2 (t - k Tp - x^2/(c H) - y^2/(c H))^2 +  2 Q (v k Tp)/(H d) + \[Tau]^2 ((v k Tp)/(H d))^2)],
     {x, \-Infinity, Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
      AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"],
  {k, -255, 255}];

ListLogPlot[Transpose[{t, p}], Joined -> True, PlotRange -> All]

ListPlot[Transpose[{t, p}], Joined -> True, PlotRange -> All]

(P0*((4*\[Pi] (\[Pi]* d^2)/^4)/\[Lambda]^2)^2*\[Lambda]^2*f*\[Tau]*KIP*\[Eta]*kf^2)/(64*\\[Pi]^2*H^2)
     Exp[-((16*\[Pi]^2*s2^2)/\[Lambda]^2)](d^4 f kf^2 KIP P0 \[Pi]^2 \[Eta] \[Tau])/
     (4 H^2 Null^8 \[Lambda]^2)E^(-((16 \[Pi]^2 s2^2)/\[Lambda]^2))(P0*((4*\[Pi] (\[Pi]* d^2)/^4)/
     \[Lambda]^2)^2*\[Lambda]^2*f*\[Tau]*KIP*\[Eta]*kf^2*c)/
    (\32*\[Pi]^2*H^3*s1^2)(c d^4 f kf^2 KIP P0 \[Pi]^2 \[Eta] \[Tau])/(2 H^3 Null^8 s1^2 \\[Lambda]^2)
z = Sum[NIntegrate[    Exp[-(5.55* r^2)/(H^2 *\[CapitalTheta]0^2) - (11.1*  r*\[CapitalTheta]s^2*Cos[\[Phi]])/
    (H *\[CapitalTheta]0^2) -       r^2/(H^2 *s1^2) - \[Pi]*       f^2*(t - k Tp - r^2/(c H) - \[Eta])^2],
    {\[Eta], -Infinity,      Infinity}, {r, 0, Infinity}, {\[Phi], 0, 2*\[Pi]},   MaxRecursion -> 40,
     AccuracyGoal -> 60,     Method -> "AdaptiveMonteCarlo"],
  {k, -255, 255}];
ListPlot[Transpose[{t, z}], Joined -> True, PlotRange -> All]

Questions concerning continuous time:
1. Why doesn't the graph show anything?
 http://zalil.ru/34911400/5efb62f5.52dc96a8/7.nb

7.nb
 c = 3*10^8;
 H = 7.2*10^5;
 \[Tau] = 47*10^-6;
 Fd = 7200;
 Tp = 50*10^-6;
 f = 3.5*10^8;
 v = 7500;
 Q = 18000;
 L = 100;
s2 = 0.003;
s1 = 2;
\[CapitalTheta]0 = 0.018;
m = 256;
\[Lambda] = 0.022;
n = 3;
lm = 25;
d = 1.2;
\[Lambda] = 0.022;
n = 3;
p[t] := Sum[
   NIntegrate[
    Exp[(\[Pi]*Fd^2*Tp^2*(m - Abs[k])^2*x^2)/
       H^2 - (n^2*\[Pi]^2*d^2*y^2)/(\[CapitalTheta]0^2*\[Lambda]^2*
         H^2*141^2) - (100*lm^2*x^2*s1^2)/(H^2*L^2*\[Pi]^2*
         s2^2) - (100*lm^2*y^2*s1^2)/(H^2*L^2*\[Pi]^2*
         s2^2) - (5.55*(y^2 +
           x^2))/(H^2*\[CapitalTheta]0^2) - \[Pi]*(f^2*(tt - k*Tp -
             x^2/(c*H) - y^2/(c*H))^2 +
         2*Q*(v*k*Tp)/(H*
             d) + \[Tau]^2*((v*k*Tp)/(H*d))^2)], {x, -Infinity,
     Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
    AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"], {k, -255,
    255}];
Plot[p[t_], {t, -10^(-8), 10^(-6)}, PlotPoints -> 30,
MaxRecursion -> 0]


I would like to correct it with your help.
POSTED BY: salangnew
6 Replies
Sort By:
salangnew
salangnew
Posted 11 years ago
I corrected. Can you put somewhere your variant of file?
 7. nb
 
 c = 3*10^8;
 H = 7.2*10^5;
 \[Tau] = 47*10^-6;
 Fd = 7200;
 Tp = 50*10^-6;
 f = 3.5*10^8;
 v = 7500;
Q = 18000;
L = 100;
s2 = 0.003;
s1 = 2;
\[CapitalTheta]0 = 0.018;
m = 256;
\[Lambda] = 0.022;
n = 3;
lm = 25;
d = 1.2;
\[Lambda] = 0.022;
n = 3;
p[t] := Sum[
   NIntegrate[
    Exp[(\[Pi]*Fd^2*Tp^2*(m - Abs[k])^2*x^2)/
       H^2 - (n^2*\[Pi]^2*d^2*y^2)/(\[CapitalTheta]0^2*\[Lambda]^2*
         H^2*141^2) - (100*lm^2*x^2*s1^2)/(H^2*L^2*\[Pi]^2*
         s2^2) - (100*lm^2*y^2*s1^2)/(H^2*L^2*\[Pi]^2*
         s2^2) - (5.55*(y^2 +
           x^2))/(H^2*\[CapitalTheta]0^2) - \[Pi]*(f^2*(t - k*Tp -
             x^2/(c*H) - y^2/(c*H))^2 +
         2*Q*(v*k*Tp)/(H*
             d) + \[Tau]^2*((v*k*Tp)/(H*d))^2)], {x, -Infinity,
     Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
   Method -> "AdaptiveMonteCarlo"], {k, -255,
    255}];
Plot[p[t_], {t, -10^(-8), 10^(-6)}, PlotPoints -> 30,
MaxRecursion -> 0]
 9.nb
 c = 3*10^8;
 H = 7.2*10^5;
 \[Tau] = 47*10^-6;
 Fd = 7200;
 Tp = 50*10^-6;
 f = 3.5*10^8;
 v = 7500;
 Q = 18000;
L = 100;
s2 = 0.003;
s1 = 2;
\[CapitalTheta]0 = 0.018;
m = 256;
\[Lambda] = 0.022;
n = 3;
lm = 25;
d = 1.2;
\[Lambda] = 0.022;
n = 3;
P0 = 25;
\[Eta] = 0.8;
KIP = 0.4;
kf = Sqrt[0.5];
\[CapitalTheta]s = 1;
t = Sort@{-10^-8, -20^-9, -10^-9, 0, 10^-11, 10^-10, 10^-9, 10^-8,
    20*10^-7, 2*10^-6, 10^-6};

p = Sum[NIntegrate[
    Exp[(\[Pi] Fd^2 Tp^2 (m - Abs[k])^2 x^2)/
       H^2 - (n^2 \[Pi]^2 d^2 y^2)/(\[CapitalTheta]0^2 \[Lambda]^2 \
H^2*141^2) - (100 lm^2 x^2 s1^2)/(H^2 L^2 \[Pi]^2 s2^2) - (100 lm^2 \
y^2*s1^2)/(H^2 L^2 \[Pi]^2 s2^2) - (5.55 (y^2 +
           x^2))/(H^2 \[CapitalTheta]0^2) - \[Pi] (f^2 (t - k Tp -
             x^2/(c H) - y^2/(c H))^2 +
         2 Q (v k Tp)/(H d) + \[Tau]^2 ((v k Tp)/(H d))^2)], {x, \
-Infinity, Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
    Method -> "AdaptiveMonteCarlo"], {k, -255,
    255}];

ListLogPlot[Transpose[{t, p}], Joined -> True, PlotRange -> All]

(\[Pi]*d^2)/^4
POSTED BY: salangnew
The 9.nb code shown has two places with Null^8. 
You should check the original notebook for extra spaces or something needing retyping.
POSTED BY: Bruce Miller
Bill Simpson
Bill Simpson
Posted 11 years ago
I apologize for the software used to post and display messages here. Each time I try to create more than one box containing code and put text before and between the boxes to explain the code I fail. Thus I posted a message trying to communicate to you in spite of these problems.

I tried to explain to you to look at the data first, before typing Plot and wondering why you cannot see your points.

and then

Plot the data with larger points so that you could see where they are.

All of that was dealing with 9.nb only.

Again I apologize for not being able to be clear.

Now in your 7.nb notebook  you have a variable tt which has not been given a value. NIntegrate demands that all variables except the variable of integration be assigned constant numeric values. Without assigning a numeric value to tt this will fail.

Some of your coefficients in 7.nb have decimal points. These tell Mathematica that the precision of those coefficients is defined by how many digits follow that decimal point. Thus you have coefficients with only one or two digits of precision, but you ask NIntegrate to find a solution with 60 digits of precision. This will likely fail. If you provide all your coefficients with 60 or more digits of precision this will have more of a chance of success, but success is not certain.

Now to try to address your question 3 and 4.

3 is perhaps due to lack of precision rounding off to zero. But even when I replace all your low precision approximate numbers with exact rationals I fail to get 60 digit precision results.

4 NIntegrate is performing an approximate decimal calculation. Existance of an analytic function may be a completely different issue.

NOTE: ALL this work was done by downloading your notebooks from your host yesterday and and again downloading your revised notebooks today. Mathematica code pasted into the text fields here for posting, perhaps unless specifically turned into InputForm, will often be corrupted by the posting process and not be able to dependably be scraped-n-pasted back into Mathematica without errors or meaningless results.

I don't think I can overcome your precision issues. Perhaps someone else can find a solution.
POSTED BY: Bill Simpson
salangnew
salangnew
Posted 11 years ago
I didn`t understand what you mean (* and then *) ?
The calculation of 7.nb doesn`t work even with definition  Plot[ p].
What`s about these 2 questions:
3.  while calculating coefficient  (at the bottom of file) the program says about dixision by 0. But it is not possible, as in Mathcad I get the result of about 10^-12
4. function z is not calculated, but it exists analytically ?
POSTED BY: salangnew
Bill Simpson
Bill Simpson
Posted 11 years ago
When I get what appears to be an empty plot I look at a table of the data if it is discrete or I create a table of the data if it is continuous.

For your 9.nb when your second ListPlot is empty I look at the data to see where the points will be. Then I adjust ListPlot to make this obvious.
Transpose[{t,p}]

(* and then *)

ListPlot[Transpose[{t, p}], Joined -> True, PlotRange -> All, PlotMarkers -> {Automatic, Medium}]

For your 7.nb I am puzzled until I notice your definition p[ t]:=... and then later Plot[ p[ t_]... and realize you have a small misunderstanding about underscore.
See if you can read this and discover the problem and the solution by yourself.
http://reference.wolfram.com/mathematica/tutorial/DefiningFunctions.html 
After you have resolved that problem there will be another problem for you to deal with, but getting the first problem resolved is an essential first step.
POSTED BY: Bill Simpson
salangnew
salangnew
Posted 11 years ago
I tried to make numerical integration of function for continuous (7.nb) and discrete times (9.nb). I have the following questions concerning discrete  time:
1.  with  logarithmic scale values of time more than 4?? are not shown, but I need to get less than 1???.
2. with linear scale,  on the contrary there appears ???????? 2us, although in the original array the maximum value is 1us
3.  while calculating coefficient  (at the bottom of file) the program says about dixision by 0. But it is not possible, as in Mathcad I get the result of about 10^-12
4. function z is not calculated, but it exists analytically

http://zalil.ru/34911397/5efb62f5.52dc96a8/9.nb
9.nb
c = 3*10^8;H = 7.2*10^5;\ = 47*10^-6;Fd = 7200;Tp = 50*10^-6;f = 3.5*10^8;v = 7500;Q = 18000;L = 100;s2 = 0.003;s1 = 2;\0 = 0.018;m = 256;\ = 0.022;n = 3;lm = 25;d = 1.2;\ = 0.022;n = 3;P0 = 25;\ = 0.8;KIP = 0.4;kf = Sqrt[0.5];\s = 1;t = Sort@{-10^-8, -20^-9, -10^-9, 0, 10^-11, 10^-10, 10^-9, 10^-8,     20*10^-7, 2*10^-6, 10^-6};p = Sum[NIntegrate[    Exp[(\ Fd^2 Tp^2 (m - Abs)^2 x^2)/       H^2 - (n^2 \^2 d^2 y^2)/(\0^2 \^2 \H^2*141^2) - (100 lm^2 x^2 s1^2)/(H^2 L^2 \^2 s2^2) - (100 lm^2 \y^2*s1^2)/(H^2 L^2 \^2 s2^2) - (5.55 (y^2 +            x^2))/(H^2 \0^2) - \ (f^2 (t - k Tp -              x^2/(c H) - y^2/(c H))^2 +          2 Q (v k Tp)/(H d) + \^2 ((v k Tp)/(H d))^2)], {x, \-Infinity, Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,     AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"], {k, -255,     255}];ListLogPlot[Transpose[{t, p}], Joined -> True, PlotRange -> All]
ListPlot[Transpose[{t, p}], Joined -> True, PlotRange -> All](P0*((4*\ (\*    d^2)/^4)/\^2)^2*\^2*f*\*KIP*\*kf^2)/(64*\\^2*H^2)Exp[-((16*\^2*s2^2)/\^2)](d^4 f kf^2 KIP P0 \^2 \ \)/(4 H^2 Null^8 \^2)E^(-((16 \^2 s2^2)/\^2))(P0*((4*\ (\*    d^2)/^4)/\^2)^2*\^2*f*\*KIP*\*kf^2*c)/(\32*\^2*H^3*s1^2)(c d^4 f kf^2 KIP P0 \^2 \ \)/(2 H^3 Null^8 s1^2 \\^2)z = Sum[NIntegrate[    Exp[-(5.55* r^2)/(H^2 *\0^2) - (11.1*         r*\s^2*Cos[\])/(H *\0^2) -       r^2/(H^2 *s1^2) - \*       f^2*(t - k Tp - r^2/(c H) - \)^2], {\, -Infinity,      Infinity}, {r, 0, Infinity}, {\, 0, 2*\},     MaxRecursion -> 40, AccuracyGoal -> 60,     Method -> "AdaptiveMonteCarlo"], {k, -255, 255}];ListPlot[Transpose[{t, z}], Joined -> True, PlotRange -> All]

Questions concerning continuous time:
1. Why doesn't the graph show anything?
 http://zalil.ru/34911400/5efb62f5.52dc96a8/7.nb

7.nb
c = 3*10^8;
H = 7.2*10^5;
\ = 47*10^-6;
Fd = 7200;
Tp = 50*10^-6;
f = 3.5*10^8;
v = 7500;
Q = 18000;
L = 100;
s2 = 0.003;
s1 = 2;
\0 = 0.018;
m = 256;
\ = 0.022;
n = 3;
lm = 25;
d = 1.2;
\ = 0.022;
n = 3;
p := Sum[
   NIntegrate[
    Exp[(\*Fd^2*Tp^2*(m - Abs)^2*x^2)/
       H^2 - (n^2*\^2*d^2*y^2)/(\0^2*\^2*
         H^2*141^2) - (100*lm^2*x^2*s1^2)/(H^2*L^2*\^2*
         s2^2) - (100*lm^2*y^2*s1^2)/(H^2*L^2*\^2*
         s2^2) - (5.55*(y^2 +
           x^2))/(H^2*\0^2) - \*(f^2*(tt - k*Tp -
             x^2/(c*H) - y^2/(c*H))^2 +
         2*Q*(v*k*Tp)/(H*
             d) + \^2*((v*k*Tp)/(H*d))^2)], {x, -Infinity,
     Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
    AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"], {k, -255,
    255}];
Plot[p[t_], {t, -10^(-8), 10^(-6)}, PlotPoints -> 30,
MaxRecursion -> 0]


I would like to correct it with your help.
POSTED BY: salangnew
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract