0
$\begingroup$

I need to find an integral, in which there is a function

YO [x_] = 
  Sum[E^(((I h PI)/(λ f))x^2) Subscript[Y, h][x, z = 13.5], {h, - 5, - 1}] /. s

where Subscript[Y, h][x, z] is a solution of partial differential equation.

I run YO [3] in the positive interval of {x, 0, 30} and find that I can get a numerical value. After that, I want to extend the domain of YO to {x, -30, 30}. The value of negative interval upper function and positive interval are symmetric. So I manually defined a piecewise function, and tried it with Piecewise and If functions. When the function YO[x] is in the negative interval, the value is equal to YO[Abs[x]], but I found that after calculation, running YO[3] can not get the numerical value, so I can only return YO[3] itself.

Based on the above situation, it is OK for me to integrate with the non expanded YO[x], but if I integrate with the expanded YO[x] with Piecewise or If function, I get the error:

the sampling point is non numerical at a certain point.

At the same time, It also report errors when I find the maximum value of the integral result g[xo, zo] by FindArgMax, Whether it is the function of YO before expansion or the function of YO after integration, It will report an error that the function value is not a real number. But my function value is the square after taking the modulus, which must be a real number. The relevant knowledge that I have been looking for for a long time these days has not been solved, so I would like to know what's wrong with my code.

In addition, I can get the specific value by directly integrating the function YO[x] which is not expanded, but I can't get the numerical value by integrating the function YO[x] which is expanded by Piecewise and If, so I can only get the function formula which contains subscript [Y, h][x,z].

This is my program for solving partial differential equations, and it is proved that the solution is completely correct.

ClearAll["Global`*"];
θ = 0;
λ = 1.2398/(19.5 10^3);
f = 4.72 10^3;
Subscript[δ, 1] = 1.274/10^6;
Subscript[δ, 2] = 4.304/10^6;
Subscript[bt, 1] = 5.254/10^9;
Subscript[bt, 2] = 2.435/10^7;
χ1 = -2 Subscript[δ, 1] + 2 I Subscript[bt, 1];
χ2 = -2 Subscript[δ, 2] + 2 I Subscript[bt, 2];
Δχ = χ1 - χ2;
k = 2 (π/λ);
Table[β[b] = -((2 (b x) Sin[θ])/f) - ((b x)/f)^2, {b, -5, 5}];
Table[χ[a] = (Δχ (1 - (-1)^Abs[a]))/(2 I a π), {a, -10, -1}];
Table[χ[a] = (Δχ (1 - (-1)^Abs[a]))/(2 I a π), {a, 1, 10}];
Table[χ[a] = (χ1 + χ2)/2, {a, 0, 0}];

eqns = 
  Join[
    Table[
      (Sin[θ] + (h x)/f) D[Subscript[Y, h][x, z], x] + 
         Cos[θ] D[Subscript[Y, h][x, z], z] == 
      ((I Pi)
        (β[h] Subscript[Y, h][x, z] +
           Sum[χ[h - l] Subscript[Y, l][x, z], {l, -5, 5}]))/λ, 
      {h, -5, 5}], 
    Table[Subscript[Y, h][x, 0] == If[h == 0, 1, 0], {h, -5, 5}]]; 

s = 
  NDSolve[eqns,
    Table[Subscript[Y, h][x, z], {h, -5, 5}],
    {x, 0, 30}, {z, 0, 30},
    Method ->
      {"MethodOfLines",
       "SpatialDiscretization" ->
         {"TensorProductGrid", 
          "MaxPoints" -> 100, "MinPoints" -> 100, "DifferenceOrder" -> 4},
      "TemporalVariable" -> z}];

This is my program to define YO[x]

YO[x_?NumericQ] :=
  Piecewise[
    {{Sum[
        E^((I h π)/(λ f) x^2) Subscript[Y, h][x, z = 13.5],
        {h, -5, -1}],
      20 >= x >= 0}, 
     {Sum[
        E^((I h π)/(λ f) x^2) Subscript[Y, h][(Abs[x]), z = 13.5], 
        {h, -5, -1}],
      0 > x >= -20}},
    {x, -20, 20}] /. s

And this is what I want to do: integrate, find the maximum of g[xo, zo] and plot the density around zo at the maximum

f1[xo_?NumericQ, zo_?NumericQ] := 
  NIntegrate[(E^(I k zo) YO[x] E^(((I k) x^2)/(2 zo))
        E^(((I k) xo^2)/(2 zo))
        E^(-(((I k) x xo)/zo)))/(I λ zo^0.5), {x, -20, 
     20}] /. s;

g[xo_?NumericQ, zo_?NumericQ] := Abs[f1[xo, zo]^2];

FindArgMax[
  {g[xo, zo], -20 <= xo <= 20 && 3000 <= zo <= 5000}, 
  {{xo, 0}, {zo, 4000}}]

DensityPlot[g[xo, zo], {xo, -20, 20}, {zo, 3000, 5000}, 
  PlotLegends -> Automatic]

But in fact, both the NIntegrate and FindArgMax will give me an error. The message from NIntegrate tells me that some values are not numerical, and the message from FindArgMax tells me that some values are not real. But it's obvious that the absolute value followed by the square must be real and numerical.

In addition, when I integrate with the unexpanded YO[x], I can integrate and draw a graph, but I never succeed in finding the maximum value.

This is my unexpanded YO[x], which is defined over {x, 0, 30}. With this function, I can run YO[3] to calculate the specific value, but if I expand the domain to {x, - 30,30}, which is the piecewise function in the previous code, I can only get itself by running YO[3].

YO[x_]=
  Sum[E^(((I h Pi)x^2)/(λ f)) Subscript[Y, h][x, z = 13.5], 
  {h, -5, -1}] /. s;

I wonder if this is the reason why plotting after integration produces complaints about non-numerical values, but no matter which YO[x] I use to calculate the maximum of g[xo,zo], it always complains about non-numerical values.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Try the following toy example and think about what's wrong: s = {{aaa -> 1}}; (aaa /. s) > 2 $\endgroup$
    – xzczd
    Jan 8, 2021 at 12:17
  • $\begingroup$ Thank you for your reply. So sometimes the output result is {a}, which will be considered as non-numerical in some functions, because it is in the form of a list? $\endgroup$
    – shrocat
    Jan 9, 2021 at 8:40
  • $\begingroup$ Yes. Another example: FindArgMax[{{-2 x^2 - 3 x + 5}}, x]. $\endgroup$
    – xzczd
    Jan 9, 2021 at 9:18
  • $\begingroup$ I see, Findargmax can calculate the list, that is {eq}, but {{}} is not a list or a number, so it cannot be calculated. @xzczd $\endgroup$
    – shrocat
    Jan 9, 2021 at 9:27
  • $\begingroup$ I don't think so. FindArgMax[{-2 x^2 - 3 x + 5}, x] works probably because it's treated as a special case of the last documented syntax FindArgMax[{f,cons},{x,y,…}]. $\endgroup$
    – xzczd
    Jan 9, 2021 at 9:37

1 Answer 1

Reset to default
0
$\begingroup$

Update

Your piecewise definition for YO is improperly written. However, I would recommend not using Piecewise at all. I would rewrite your code like this:

θ = 0;
λ = 1.2398/(19.5 10^3);
f = 4.72 10^3;
Subscript[δ, 1] = 1.274/10^6;
Subscript[δ, 2] = 4.304/10^6;
Subscript[bt, 1] = 5.254/10^9;
Subscript[bt, 2] = 2.435/10^7;
χ1 = -2 Subscript[δ, 1] + 2 I Subscript[bt, 1];
χ2 = -2 Subscript[δ, 2] + 2 I Subscript[bt, 2];
Δχ = χ1 - χ2;
k = 2 (π/λ);
Table[β[b] = -((2 (b x) Sin[θ])/f) - ((b x)/f)^2, {b, -5, 5}];
Table[χ[a] = (Δχ (1 - (-1)^Abs[a]))/(2 I a π), {a, -10, -1}];
Table[χ[a] = (Δχ (1 - (-1)^Abs[a]))/(2 I a π), {a, 1, 10}];
Table[χ[a] = (χ1 + χ2)/2, {a, 0, 0}];

eqns = 
  Join[
    Table[
      (Sin[θ] + (h x)/f) D[Subscript[Y, h][x, z], x] + 
       Cos[θ] D[Subscript[Y, h][x, z], z] == 
      ((I Pi) 
        (β[h] Subscript[Y, h][x, z] + 
          Sum[χ[h - l] Subscript[Y, l][x, z], {l, -5, 5}]))/λ,
      {h, -5, 5}],
    Table[Subscript[Y, h][x, 0] == If[h == 0, 1, 0], {h, -5, 5}]];

The following contains modifications to your code.

s =
  NDSolve[eqns,
    Table[Subscript[Y, h], {h, -5, 5}], {x, 0, 30}, {z, 0, 30},
    Method ->
      {"MethodOfLines",
       "SpatialDiscretization" ->
         {"TensorProductGrid", "MaxPoints" -> 100,
          "MinPoints" -> 100, "DifferenceOrder" -> 4},
       "TemporalVariable" -> z}];

Clear[YO]
YO[x_?NumericQ /; (0 > x >= -20)] := YO[-x]
YO[x_?NumericQ] = 
  First[Sum[E^((I h π)/(λ f) x^2) Subscript[Y, h][x, 13.5], {h, -5, -1}] /. s];

Then YO behaves as a even function over the domain [-20, 20]:

YO /@ {-20, -3, 3, 20} // Column

result

I think there are (possibly many) other problems with your code, but maybe this will help you to move on and correct them.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Thank you for your answer, which is of great help to me. You use the First to change the output {a} into a, Let me know what I made by outputting errors. Although YO[x] has been able to get the value, but I still have some things I don't understand. For example, in NDSolve, you replace subscript [Y, h] [x, z] with subscript [Y, h], and then you can get the value of YO [x]. I still don't understand it after I looked up some information. $\endgroup$
    – shrocat
    Jan 9, 2021 at 5:51
  • $\begingroup$ @shrocat Try the following toy example and think about what's wrong: func[x_]=x+1; f[1] /. f[x]->func[x], then try the following and think about why it works: func[x_]=x+1; f[1] /. f -> func . $\endgroup$
    – xzczd
    Jan 9, 2021 at 7:31
  • $\begingroup$ Hello,ba zhu. I feel that I understand something, but I don't seem to understand it. The following function /. Replaces the head, so the numerical result can be obtained. If the previous replacement is replaced by func[x_ ] = x + 1; f[1] /. f[x_ ]- > func [x] or func[x] = x + 1; f [1] /. F [x]_ ]- > func [x] can also get numerical results, so I think it should be a variable problem? $\endgroup$
    – shrocat
    Jan 9, 2021 at 8:38
  • $\begingroup$ @shrocat You need to add @xzczd in the comment or I won't get the reminder. Then, something is missing in your second sample, please double check it. $\endgroup$
    – xzczd
    Jan 9, 2021 at 9:24
  • $\begingroup$ @xzczd Oh,I made some mistakes.It should be func[x] = x + 1; f[1] /. f[x_] -> func[x] and func[x_] = x + 1; f[1] /. f -> func. $\endgroup$
    – shrocat
    Jan 9, 2021 at 9:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.