What is the correct ContourPlot here? Where should I assign my parameters?

GROUPS:
 I am currently trying to make a contour plot of a transmission function. It is a very naughty function as some parts will tend to infinity and others to zero, however as a whole it will (in the allowed cases at least) stay between 0 and 1, as it should. This function depends on five parameters and these parameters are used to define expressions to simplify the function because it is huge. Just as an example, this is what my code would look like: V = 5; U = 4; h = 1.6; v = 0.6; (*assigning a numerical value to the parameters*) K = U + V +hv(a + b) (*defining an expression where 'a' and 'b' are variables.*) T[a_, b_] = K + hv (*defining a function*) ContourPlot[T[a,b], {a, 0.0, 1.0}, {b, 0.0, 1.0}] Now the problem comes when I change the order in which I assign values to my parameters. If I do it at the beginning of my program I'll have one ContourPlot (1), however if I assign values to my parameters just before I order the plot; for example: K = U + V +hv(a + b) T[a_, b_] = K + hv V = 5; U = 4; h = 1.6; v = 0.6; ContourPlot[T[a,b], {a, 0.0, 1.0}, {b, 0.0, 1.0}] I'll get something very different (2).So I tried simplifying this last expression before the plot and after assigning the values: K = U + V +hv(a + b) T[a_] = K + hv V = 5; U = 4; h = 1.6; v = 0.6; T[a_, b_] = Simplify[ T[a_, b_] ] ContourPlot[T[a,b], {a, 0.0, 1.0}, {b, 0.0, 1.0}] and the result was again (1).(1) (2)I read the Evaluation procedure, and it seems logical that Mathematica was simplifying the expressions as much as it could with the parameters given. When the parameters were not given numerical values, Mathematica did not simplify and ended up with a more complex function and given its naughtiness, I assume some information may change in the process.The question is, what is the correct plot/process? How can identify the correct procedure in these cases? Thank you in advance.Edit: I posted two notebooks where this same thing happens. Attachments:
1 month ago
23 Replies
 I think you'd get more and better responses if you supplied a complete working example. And admittedly I can be a bit dense as I'm not at all certain what the issue is that you're having trouble with.
1 month ago
 Thanks for the comment Jim. Now I can't put the code but I'll try to reproduce this problem with a simpler function and I'll post it as soon as I can.
1 month ago
 Jim Baldwin 1 Vote I still don't understand what the issue is but a wild guess is that you don't like the fact that parts of one of the contour plots is "white". Usually to fix that one can use PlotRange -> All.
1 month ago
 I'll admit I hate the white spots, haha, but that's not the problem. The problem is that both contour plots that appear above come from the same function (and the same ContourPlot instructions). The order in which I assign values to my parameters alters the plot.
1 month ago
 Marco Thiel 1 Vote Yes, and if that is the problem also ClippingStyle -> False should work: Row[{ContourPlot[1/(Sin[x^2 + y^2]^2 + 0.02), {x, -2, 2}, {y, -2, 2}, ImageSize -> Medium], ContourPlot[1/(Sin[x^2 + y^2]^2 + 0.02), {x, -2, 2}, {y, -2, 2}, PlotRange -> All, ImageSize -> Medium], ContourPlot[1/(Sin[x^2 + y^2]^2 + 0.02), {x, -2, 2}, {y, -2, 2}, ClippingStyle -> False, ImageSize -> Medium]}] But I somehow feel that that is not the actual problem. Cheers,Marco
1 month ago
 Gianluca Gorni 1 Vote Simplify[T[a_]] does not simplify the definition of T. For that you must insert Simplify into the definition: T[a_]=Simplify[k+hv].
1 month ago
 Thank you for the remark Gianluca, I have already edited that mistake.
1 month ago
 Hans Dolhaine 1 Vote And what is hv?
1 month ago
 Planck's constant and escape velocity but I'll edit that to make it clearer, thanks for the reply :)
1 month ago
 Hans Dolhaine 1 Vote I am a bit puzzled. You define T[a_] = ..... and you do a ContourPlot[ T[a, b].Your definition of T is not a SetDelayed : T[a_] := expr and T[a,b] is nowhere defined. So the ContourPlot shouldN#t work at all.Furthermore your function is linear in a and b. When I put things together I arrive at Clear[TT] TT[a_, b_] := U + V + h v (a + b) + h v and this shouldn't give the ContourPlots you show us. Perhaps it is a good idea when you show a notebook with a bit more code.And perhaps is your problem that T doesn't know the right parameters. Perhaps you should define Clear[TT] TT[a_, b_, U_, V_, h_, v_] := U + V + h v (a + b) + h v but this is still linear in a and b
1 month ago
 Marco Thiel 1 Vote This time I was first (even though my post appears below)!;-)Marco
1 month ago
 Hans Dolhaine 1 Vote smile
1 month ago
 Marco Thiel 1 Vote The thing that intrigues me most is what exactly the definition T[a_] = K + hv has got to do with the plot - it is a function of one variable. There ContourPlot[T[a,b], {a, 0.0, 1.0}, {b, 0.0, 1.0}] T is a function of two variables. So the definition is inconsequential and should give an empty image, should it not? It therefore appears that T as a function of two variables must be defined somewhere else in the code/notebook. It would be useful to have that definition.Cheers,Marco
1 month ago
 Yes, both Hans and Marco are totally right. I made a mistake in the example. I apologise for the confusion. I uploaded a 2 working notebooks, however, so it becomes more clear. I also edited that mistake now. Thanks again,Edit: The notebooks contain a plot of what should be the same function. They vary only on the place where I assign the values of my parameters.
1 month ago
 Marco Thiel 1 Vote The notebook ContourPlot 2.nbgives the same as the first if you comment out the first three lines (the Clear commands) after second evaluation.ClearAll doesn't appear to do anything anyway.I think.Cheers, Marco
 Marco Thiel 1 Vote If I run (*ClearAll["Global*"]*) (*ClearAll["Global*"] ClearAll[Evaluate[$Context<>"*"]]*) ClearAll["Global*"]; energy = 0.027*eV; eV = 1/27.2113845; nm = 1/0.052917721092; a0 = 4.65; n = 0; M = 4; W = a0*(3*M + 1); e = 1; h = 1; v = 0.003*137.036; kmodC = (energy - e*VC)/(h*v); kmodRG = (energy - e*VRG)/(h*v); kmodRB = (energy - e*VRB)/(h*v); kC = (kmodC^2 - kn^2)^0.5; kRG = (kmodRG^2 - kn^2)^0.5; kRB = (kmodRB^2 - kn^2)^0.5; zC = (kn + I*kC)/((kn^2 + kC^2)^0.5); zRG = (kn + I*kRG)/((kn^2 + kRG^2)^0.5); zRB = (kn + I*kRB)/((kn^2 + kRB^2)^0.5); T22 = ( E^(-2 I a (kC + kRG) + I b (2 kRG + kRB)) ((zC zRG - 1)^2 E^( 4 I (a - b) kRG) (-1 + zRB zRG)^2 + 2 E^(2 I (a - b) kRG) (-zC (1 + zRG^2) + zRG + zC^2 zRG)))/((-1 + zRG^2)^2); T[kn_, a_] = 1/(T22*(T22)\[Conjugate]); (*Parameters*) VRG = 0.100*eV; VRB = 0.001*eV; VC = 0.001*eV; b = (a + d); d = 50*nm; (*Parameters*) T22post = ( E^(-2 I a (kC + kRG) + I b (2 kRG + kRB)) ((zC zRG - 1)^2 E^( 4 I (a - b) kRG) (-1 + zRB zRG)^2 + 2 E^(2 I (a - b) kRG) (-zC (1 + zRG^2) + zRG + zC^2 zRG)))/((-1 + zRG^2)^2); Tpost[kn_, a_] = 1/(T22*(T22)\[Conjugate]); Evaluate[T[kn*(1/nm), (a/2)*nm] === Tpost[kn*(1/nm), (a/2)*nm]] It turns out that the T and Tpost are not identical as Hans says. The first bits of the functions are a bit different:TPostat the very end the "a" has a different factor.In principle it is now easy to see what is happening. Just compare TreeForm[Tpost[kn*(1/nm), (a/2)*nm]] to TreeForm[T[kn*(1/nm), (a/2)*nm]] I agree that the delayed definition would probably be more appropriate.Best wishes,Marco Answer 1 month ago  Ok, so what I understand is that, there is something odd happening, the function is definitely not the same depending on how I assign values to my parameters as Marco clearly demonstrated. Both Marco and Hans, you are advicing me to hold any evaluation until the plot. Did I get it right?You are also suggesting this procedure would give the correct answer. In your opinion, can I really trust this is the case or should I use a different program just to be sure, like Matlab; or will I find the same issue there?Thanks. Answer 1 month ago  Hans Dolhaine 1 Vote I think the problem is solved. It seems that there were problems with handing the parameters to the function(s). Look at my code ( a modificaton of Marco's code). It gives the same functions regardless where you do the assignments. And in my opinion Mathematica will do the job. Answer 1 month ago  Ok, then I'll rest assured that delaying is the correct procedure, Hans. Thank you very much for your time. :) Answer 1 month ago  Thank your very much Marc. I'll rest assured that delaying the definition will return the correct plot. Answer 1 month ago  Hans Dolhaine 1 Vote I had another look at the problem. With the code below (copy of Marco's code above with slight modifications) everything seems to be ok. I have the idea that during the definitions of the functions Mathematica does some evaluations, and the result of these depend on whether numbers are already assigned or not. But this is just an assumption. Note that I use Setdelayed and plug in the parameters in T22 (and T22post) as variables of these functions. The code give the same results for T and Tpost even if you switch all the Clear-statements at the beginning on.@Marco: Which is the reason you don't use T22post inTpost[kn, a] = 1/(T22*(T22)[Conjugate]);And here my code: In[1]:= (*ClearAll["Global*"]*)(*ClearAll["Global*"] \ ClearAll[Evaluate[$Context<>"*"]]*)ClearAll["Global*"]; energy = 0.027*eV; eV = 1/27.2113845; nm = 1/0.052917721092; a0 = 4.65; n = 0; M = 4; W = a0*(3*M + 1); e = 1; h = 1; v = 0.003*137.036; kmodC = (energy - e*VC)/(h*v); kmodRG = (energy - e*VRG)/(h*v); kmodRB = (energy - e*VRB)/(h*v); kC = (kmodC^2 - kn^2)^0.5; kRG = (kmodRG^2 - kn^2)^0.5; kRB = (kmodRB^2 - kn^2)^0.5; zC = (kn + I*kC)/((kn^2 + kC^2)^0.5); zRG = (kn + I*kRG)/((kn^2 + kRG^2)^0.5); zRB = (kn + I*kRB)/((kn^2 + kRB^2)^0.5); T22[kn_, a_] := Evaluate[(E^(-2 I a (kC + kRG) + I b (2 kRG + kRB)) ((zC zRG - 1)^2 E^(4 I (a - b) kRG) (-1 + zRB zRG)^2 + 2 E^(2 I (a - b) kRG) (-zC (1 + zRG^2) + zRG + zC^2 zRG)))/((-1 + zRG^2)^2) /. b -> a + d]; T[kn_, a_] = 1/(T22[kn, a]*(T22 [kn, a])\[Conjugate]); (*Parameters*) VRG = 0.100*eV; VRB = 0.001*eV; VC = 0.001*eV; b = (a + d); d = 50*nm; (*Parameters*) T22post[kn_, a_] := Evaluate[(E^(-2 I a (kC + kRG) + I b (2 kRG + kRB)) ((zC zRG - 1)^2 E^(4 I (a - b) kRG) (-1 + zRB zRG)^2 + 2 E^(2 I (a - b) kRG) (-zC (1 + zRG^2) + zRG + zC^2 zRG)))/((-1 + zRG^2)^2) /. b -> a + d]; Tpost[kn_, a_] = 1/(T22post[kn, a]*(T22post[kn, a])\[Conjugate]); Evaluate[T[kn*(1/nm), (a/2)*nm] === Tpost[kn*(1/nm), (a/2)*nm]] Out[30]= True Example In[30]:= T[1, 1] Tpost[1, 1] Out[30]= 2.927595687384534*10^-822 + 0.*10^-850 I Out[31]= 2.92759568738453*10^-822 + 0.*10^-850 I In[32]:= T22[1, 1] T22post[1, 1] Out[32]= 5.844460807504422*10^410 - 3.400313161721839*10^397 I Out[33]= 5.84446080750442*10^410 - 3.40031316172184*10^397 I