# Strange behavior when replacing variables by numerical values [duplicate]

I have a rather complicated function with parameters {a, b, c, d, e, f, k}, and I'd like to know its behavior as a function of k alone given other parameters, so I try the following code:

myFunction = D[ComplexExpand[Arg[((a^2 - 2 (2 k + I e - 2 c) (k + d - c))^2/((a^2 - 2 (2 k + I (b + e + 2 I c)) (k + d - c))^2 + 4 E^(2 I k f) b^2 (k + d - c)^2))], TargetFunctions -> {Re, Im}] // Simplify, k];
myFunctionWithValues[k_] := Evaluate[myFunction /. {b -> 10, a -> 2, c -> 100, d -> 0, e -> 1, f -> Pi/100} // FullSimplify]
LogPlot[myFunctionWithValues[k], {k, 95, 105}, PlotRange -> All]


and the outcome looks very bizarre. It can't be correct. To comfirm this, for example, I look at the value of myFunctionWithValues at the center k=100 by 1) directly replacing k by 100 and 2) by taking the limit k->100, and they disagree with each other:

myFunctionWithValues[k] /. k -> 100 // N


315120.

Limit[myFunctionWithValues[k], k -> 100] // N


10.

Apparently the latter gives a more reasonable result (which is checked by other approaches). Moreover, if I replace the value of all parameters (including k) in the very beginning, Mathematica also yields the same result:

myFunction /. {b -> 10, a -> 2, k -> c, d -> 0, e -> 1, f -> Pi/100} // Simplify


10

myFunction /. {b -> 10, a -> 2, c -> 100, d -> 0, e -> 1, f -> Pi/100, k -> 100} // N


10.

And the correct plot can be given by the following two ways:

(Limit[myFunctionWithValues[k], k -> #]) & /@ Range[95, 105, 1/10] //N;
ListPlot[%, PlotRange -> All, Joined -> True, AxesOrigin -> {100, 0}, DataRange -> {95, 105}, PlotMarkers -> Automatic]


Plot[myFunction /. {b -> 10, a -> 2, c -> 100, d -> 0, e -> 1, f -> Pi/100}, {k, 95, 105}, PlotRange -> All, AxesOrigin -> {100, 0}]


Question: What's going on here? Why can't I simply plot myFunctionWithValues but have to take the limit of k at each point? More precisely, what's the difference between using the replacement rules and taking limits? I thought it should be the same given that in my case there's no singularity causing trouble, otherwise Mathematica would give me an error message.

-
I suspect there is some kind of singularity. Look at myFunctionWithValues[100.03] // N and then at 100.02 and at 100.01, the values jump all over the place. You are dividing by things, and probably one of these is very small. –  bill s Jul 17 '13 at 21:01
@bills, that was my first guess too. But there shouldn't be any singularity! For example, Plot[myFunction /. {b -> 10, a -> 2, c -> 100, d -> 0, e -> 1, f -> Pi/100}, {k, 95, 105}, PlotRange -> All] also gives the correct plot which doesn't blow up. So if there is any singular behavior in myFunctionWithValues, it must be an artificial one. –  Leo Fang Jul 17 '13 at 21:08
Leo, I have marked this question as a duplicate because I believe the problem is identical. Simply, by default plotting is done with machine precision which isn't sufficient in this case. If you try LogPlot[myFunctionWithValues[k], {k, 95, 105}, PlotRange -> All, WorkingPrecision -> 50] you will get the plot that you desire. –  Mr.Wizard Jul 18 '13 at 1:28
Aha! Here comes again the precision issue! So glad you pointed out the other post to me! (By the way if you were not the one who answered that post, I would doubt that how come could someone come up a link with it...It seems unlikely for me to find that post simply by searching.) Thanks @Mr.Wizard! –  Leo Fang Jul 18 '13 at 2:48
Leo, I agree it is often hard to find past questions/answers, even when I know they exist! Sometimes Google is better than the SE (integrated) search, though the SE search works pretty well these days. Also, questions that get marked as duplicates, such as this one, serve the important purpose of providing additional chances to find the question. –  Mr.Wizard Jul 19 '13 at 0:16

## marked as duplicate by Mr.Wizard♦Jul 18 '13 at 1:25

I'd define your functions with all the arguments :

myFunction[a_, b_, c_, d_, e_, f_, k_] = D[ComplexExpand[
Arg[((a^2 - 2 (2 k + I e - 2 c) (k + d - c))^2/((a^2 - 2 (2 k + I (b + e + 2 I c))
(k + d - c))^2 + 4 E^(2 I k f) b^2 (k + d - c)^2))],
TargetFunctions -> {Re, Im}] // Simplify, k];

myFunctionWithValues[k_] = myFunction[2, 10, 100, 0, 1, Pi/100, k];

Plot[{myFunctionWithValues[k], Log[myFunctionWithValues[k]]}, {k, 95, 105},
PlotRange -> All]


-
thanks for the feedback. Actually I found another way to produce the plot just before you posted this answer. (Please see my edited post.) But I am still confused why I got the wrong plot in the beginning and expecting some sort of explanation. –  Leo Fang Jul 17 '13 at 22:19