Message Boards Message Boards

0
|
7418 Views
|
23 Replies
|
14 Total Likes
View groups...
Share
Share this post:

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

Posted 7 years ago

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).

Contour plot of my function of transmission when assigning values to my parameters at the beginning of my program(1) Contour plot of my function of transmission when assigning values to my parameters after definning my function and before plotting(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:
23 Replies

Dear Hans,

@Marco: Which is the reason you don't use T22post in

Tpost[kn, a] = 1/(T22*(T22)[Conjugate]);

It is difficult to come up with a credible excuse for my stupid mistake. Sorry!

Best wishes,

Marco

POSTED BY: Marco Thiel

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 in

Tpost[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
POSTED BY: Hans Dolhaine

The notebook ContourPlot 2.nb

gives 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

POSTED BY: Marco Thiel

That is not the case on my system.

I inserted a statement to print T22 in both notebooks. They are quite different.

I suggest to define functions T22 and T with ALL parameters as arguments (I won't do it, it is late now). And I would define all functions with a SetDelayed

Attachments:
POSTED BY: Hans Dolhaine

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:

T

enter image description here

Post

enter image description here

at 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

POSTED BY: Marco Thiel

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.

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.

POSTED BY: Hans Dolhaine

Ok, then I'll rest assured that delaying is the correct procedure, Hans. Thank you very much for your time. :)

Thank your very much Marc. I'll rest assured that delaying the definition will return the correct plot.

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

POSTED BY: Marco Thiel

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.

And what is hv?

POSTED BY: Hans Dolhaine

Planck's constant and escape velocity but I'll edit that to make it clearer, thanks for the reply :)

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

POSTED BY: Hans Dolhaine

This time I was first (even though my post appears below)!

;-)

Marco

POSTED BY: Marco Thiel

smile

POSTED BY: Hans Dolhaine

Simplify[T[a_]] does not simplify the definition of T. For that you must insert Simplify into the definition: T[a_]=Simplify[k+hv].

POSTED BY: Gianluca Gorni

Thank you for the remark Gianluca, I have already edited that mistake.

Posted 7 years ago

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.

POSTED BY: Jim Baldwin

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.

Posted 7 years ago

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.

POSTED BY: Jim Baldwin

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.

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]}]

enter image description here

But I somehow feel that that is not the actual problem.

Cheers,

Marco

POSTED BY: Marco Thiel
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract