# Plotting Function in a function not working

I consistently have problems using the plot tool when I want to find how some function (like maximization) varies as I change a parameter.

In this example, I want to find how a function's imaginary part obtains its maximum as a function of another parameter. So I want to plot this maximum calcuation vs the parameter that I am varying.

Here is my attempt at the code:

    NewIm[expression_] := (num = Numerator[expression];
den = Denominator[expression];
{reNum, imNum} = ComplexExpand[ReIm[num]]; {reDen, imDen} =
ComplexExpand[ReIm[den]];
imout = (imNum reDen - reNum imDen)/(reDen^2 + imDen^2));

NewRe[expression_] := (num = Numerator[expression];
den = Denominator[expression]; {reNum, imNum} =
ComplexExpand[ReIm[num]]; {reDen, imDen} =
ComplexExpand[ReIm[den]];
reout = (reNum reDen + imNum imDen)/(reDen^2 + imDen^2));

MaxFINDER[\[Chi]_] :=
NMaximize[
Simplify[
NewIm[\[Chi]]], \[CapitalOmega]c > .001 && \[CapitalOmega]d >
0, {\[CapitalOmega]c, \[CapitalOmega]d, \[CapitalDelta]p, \
\[CapitalDelta]c, \[CapitalDelta]s, \[CapitalDelta]d, \[Phi]},
MaxIterations -> 10000];

pCs\[Chi]2fun[\[Beta]_]  := (I E^(I \[Phi]) \[CapitalOmega]c)/(
8 (-(\[Gamma]a/2) + I \[CapitalDelta]a) ((\[Beta] \[Gamma]a)/2 + (
I \[Gamma]a \[CapitalDelta]c)/2 - I \[Beta] \[CapitalDelta]p - (
I \[Gamma]a \[CapitalDelta]p)/
2 + \[CapitalDelta]c \[CapitalDelta]p - \[CapitalDelta]p^2 + \
\[CapitalOmega]c^2/4)) /. {\[Gamma]a ->
1, \[CapitalDelta]a -> \[CapitalDelta]c + \[CapitalDelta]p - \
\[CapitalDelta]s};

(*level5\[Chi]2 = -(\[ExponentialE]^(\[ImaginaryI] \[Phi])(50 \
\[CapitalOmega]c \[CapitalOmega]d)/((\[Gamma]a^2+2 \[ImaginaryI] \
\[Gamma]a (50 \[CapitalDelta]c-51 \[CapitalDelta]p)+50 (4 \
\[CapitalDelta]c \[CapitalDelta]p-4 \
\[CapitalDelta]p^2+\[CapitalOmega]c^2)) ((\[Gamma]a+2 \[ImaginaryI] (\
\[CapitalDelta]c-\[CapitalDelta]p-\[CapitalDelta]s)) (\[Gamma]a+2 \
\[ImaginaryI] (\[CapitalDelta]c+\[CapitalDelta]d-\[CapitalDelta]p-\
\[CapitalDelta]s))+\[CapitalOmega]d^2)))/.{\[Gamma]a\[Rule]1, \
dephasing\[Rule]10^-4,\[CapitalDelta]a->-1(-\[CapitalDelta]p +\
\[CapitalDelta]c-\[CapitalDelta]s+\[CapitalDelta]d)};*)
(*
MaxFINDER[psC\[Chi]2]
MaxFINDER[level5\[Chi]2]*)
MaxFINDER[pCs\[Chi]2fun[.1]]
Plot[Evaluate[MaxFINDER[pCs\[Chi]2fun[bacon]]], {bacon, .1, 2},
PlotPoints -> 15]


As you can see, if I run the fuction I created, I can find the maximum values for individual values of beta, but I'm still can't figure out how to get the plot tool to actually plot these.

The plot comes out blank, and I've tried a number of different suggestions that I've seen on previous stack-exchange questions.

Any ideas?

• Try Plot[MaxFINDER[pCs\[Chi]2fun[bacon]][[1]], {bacon, .1, 2}, PlotPoints -> 15], since NMaximize returns a list, where the first part is the max value, and the second part is the value of all the parameters. – MelaGo Jul 15 '19 at 21:34
• oh woops! thank you! – Steven Sagona Jul 15 '19 at 21:35

Your code is a bit of a mess, but here is a slightly cleaner version. First of all, user-defined functions should begin with lower case letter; second, you should localise your variables; third, if your functions depend on the denominator and numerator of your arguments, you should Rationalize those, for otherwise Numerator[0.1] == 0.1 and Denominator[0.1] == 1 instead of 1 and 10; and fourth, you should use patterns to make sure the functions are only evaluated for numeric arguments. Finally, the minimisation step is rather slow, so it is better to generate points and ListPlot them instead of just Plotting the actual function.

Clear[newIm, newRe, maxFINDER, pCsχ2fun]
newIm[expression_?NumericQ] := Module[{num = Numerator[Rationalize[expression, 0]], den = Denominator[Rationalize[expression, 0]], reNum, imNum, reDen, imDen},
{reNum, imNum} = ReIm[num]; {reDen, imDen} = ReIm[den];
(imNum reDen - reNum imDen)/(reDen^2 + imDen^2)
];
newRe[expression_?NumericQ] := Module[{num = Numerator[Rationalize[expression, 0]], den = Denominator[Rationalize[expression, 0]], reNum, imNum, reDen, imDen},
{reNum, imNum} = ReIm[num]; {reDen, imDen} = ReIm[den];
(reNum reDen + imNum imDen)/(reDen^2 + imDen^2)
];

pCsχ2fun[β_?NumericQ, Δc_, Δp_, Δs_, ϕ_, Ωc_] := With[{γa = 1, Δa = Δc + Δp - Δs}, (I E^(I ϕ) Ωc)/(8 (-(γa/2) + I Δa) ((β γa)/2 + (I γa Δc)/2 - I β Δp - (I γa Δp)/2 + Δc Δp - Δp^2 + Ωc^2/4))]

maxFINDER[bacon_?NumericQ] := NMaximize[{newIm[pCsχ2fun[bacon, Δc, Δp, Δs, ϕ, Ωc]], Ωc > .001 && Ωd > 0}, {Ωc, Ωd, Δp,  Δc, Δs, Δd, ϕ}, MaxIterations -> 10000][[1]];

ListPlot[Table[{b, maxFINDER[b]}, {b, .1, 2, .1}], Joined -> True]


On a second thought, I'm not sure why you even define newRe and newIm at all: these are identical to the built-in Re and Im! Your code can be simplified a lot:

maxFINDER[β_?NumericQ] := NMaximize[{(-Ωc (β (2-8 Δp (Δc+Δp-Δs))+4 (Δc-Δp) (Δc+2 Δp-Δs)+Ωc2) Cos[ϕ]+2 Ωc (-Δc+2 β Δc+Δp+4 β Δp+4 Δc2 Δp-4 Δp3-2 β Δs-4 Δc Δp Δs+4 Δp2 Δs+(Δc+Δp-Δs) Ωc2) Sin[ϕ])/((1+4 (Δc+Δp- Δs)2) (4 (β2+(Δc-Δp)2) (1+4 Δp2)+4 (β+2 (Δc-Δp) Δp) Ωc2+Ωc4)), Ωc > .001 && Ωd > 0}, {Ωc, Ωd, Δp, Δc, Δs, Δd, ϕ}, MaxIterations->10000][[1]]


which yields the same plot as before, with no need to define any auxiliary functions.

• Thanks for the tips. I'll try to go through the cleaner code and learn how to program more efficently with what you've wrote. – Steven Sagona Jul 15 '19 at 21:56
• @StevenSagona I'm glad I could help. Out of curiosity, why did you define your own NewRe and NewIm functions? these are identical to the built-ins Re and Im, respectively! – AccidentalFourierTransform Jul 15 '19 at 22:13
• Simplifying things with big denominators appear to have some kind of bug or problem. I asked the question and you can find it here: mathematica.stackexchange.com/questions/193509/… – Steven Sagona Jul 15 '19 at 23:14
• Going through your answer in more detail, what is the reason for using With function instead of the replace feature? Maybe its that mathematica can more easily evaluate nested expressions with these features? – Steven Sagona Jul 15 '19 at 23:20
• @StevenSagona 1) When you call Re and Im in this specific problem, it is only evaluated for actual numerical expressions, not symbolic ones, so the built-ins will not hang; they will actually be much faster than your user-defined functions. 2) There is no specific reason to use With instead of replace here; it's just more natural to me (and feels more "correct", style-wise). Also, it is easier to read (although some people might find /. easier; kind of like $\int dx\ f(x)$ vs $\int f(x)dx$). – AccidentalFourierTransform Jul 15 '19 at 23:53