Parametric Plots in Mathematica

[hufflepup]

I have 2 parameters, p and q.

rmin depends on p and q in a complex way, but for a given p and q it can be solved for numerically.

I have three functions

L(p,q,rmin)
Theta(p,q,rmin)
S(p,q,rmin)

What I want to do is to plot S against L for a set of the values p,q for which Theta is fixed. Can Mathematica do this efficiently?

My code is below:



K[r_] = 1 + b^2/r^2

h[r_] = 1 - 1/r^4

gtt[r_] = -h[r]*K[r]^(-1)*(r/R)^2*(1 + b^2*r^2)*rh^2

gyt[r_] = -2*b^2*r^2*(K[r]^(-1)*(r/R)^2)*h[r]*rh^2

gyy[r_] = K[r]^(-1)*(r/R)^2*(1 - b^2*r^2*h[r])*rh^2

gzz[r_] = (r/R)^2*rh^2

grr[r_] = (R/r)^2*h[r]^(-1)
gxx[r_] = (r/R)^2*rh^2

s = (c^2 - 1)^(1/2)

ggtt[r_] = FullSimplify[c^2*gtt[r] + s^2*gzz[r]]
ggyt[r_] = FullSimplify[c*gyt[r]]
ggzt[r_] = FullSimplify[-2*c*s*(gtt[r] + gzz[r])]
ggyz[r_] = FullSimplify[-s*gyt[r]]
ggzz[r_] = FullSimplify[c^2*gzz[r] + s^2*gtt[r]]
ggyy[r_] = FullSimplify[K[r]^(-1)*(r/R)^2*(b*R)^2*(1 - h[r])]
ggrr[r_] = FullSimplify[(R/r)^2*h[r]^(-1)]
ggxx[r_] = FullSimplify[gxx[r]]


(* Note dil[r]=e^(-2\[CapitalPhi])*)

dil[r] = K[r]

c = 1.0
b = 1.0
rh = 1
R = 1
\[CapitalLambda] = 100000


xsigma[r_] = p/q*(-gtt[r]*gzz[r] + gtx[r]^2)/(-gtt[r]*gxx[r])


rsigma[r_] =
1/q/(-gtt[r]*grr[r])^(0.5)*(-gtt[r]*gzz[r] +
gtz[r]^2)^(0.5)*((-gtt[r]*gzz[r] + gtz[r]^2) (dil[r] +
p^2/(gtt[r]*gxx[r])) - q^2)^(0.5)


rmin[p_, q_?NumericQ] :=
NSolve[(-gtt[r]*gzz[r] + gtz[r]^2) (dil[r] + p^2/(gtt[r]*gxx[r])) -
q^2 == 0, r]

Theta[rmin_, p_, q_?NumericQ] :=
NIntegrate[xsigma[r]/rsigma[r], {r, rmin, \[CapitalLambda]}]/
NIntegrate[1/rsigma[r], {r, rmin, \[CapitalLambda]}]

L[rmin_, p_, q_?NumericQ] :=
2*((NIntegrate[
1/rsigma[r], {r, rmin, \[CapitalLambda]}])^2 + (NIntegrate[
xsigma[r]/rsigma[r], {r, rmin, \[CapitalLambda]}])^2)^(0.5)


S[rmin_, p_, q_?NumericQ] :=
dil[r]/q*(gtz[r]^2 - gtt[r]*gzz[r])/rsigma[r]

ParametricPlot[{L[p_, q_], S[p_, q_]}, {p, 0, 10},
AspectRatio -> 1/GoldenRatio, Frame -> True, Axes -> False]

If I wanted to do a Log plot of this with a parametric graph would this be possible in Mathematica?

 

[Dale]

I have looked over your code and have a few suggestions. First, the symbol K is a built in symbol, I would recommend using the lower-case k instead. Second, gtx and gtz are used in other expressions but not defined. Third, your use of Set (=) and DelayedSet (:=) are pretty subtle, I think it is correct, but you should be very careful and test each part to make sure each individual part works as expected.

Now, as for your actual question. The way I would approach it is as follows:
I would introduce

qOfPTheta[p_,theta_] := NSolve[theta==Theta[p,q],q][[1,2]]

then set theta to be whatever value you are interested in and use it in

ParametricPlot[{L[p, qOfPTheta[p,theta]], S[p, qOfPTheta[p,theta]]}, {p, 0, 10},
AspectRatio -> 1/GoldenRatio, Frame -> True, Axes -> False]

 

[hufflepup]

Thanks for your reply. I have made the changes you suggested, however I am getting the following error messages when I run my new code:

Inverse functions are being used by Solve, so some solutions may not \
be found; use Reduce for complete solution information.

and

"Part 2 of \
{q->InverseFunction[0.8,2,2][0.*10^-308,0.\
8000000000000000000000000000000]} does not exist."



If I replace the NSolve's with Reduce then I get

Reduce::inex: Reduce was unable to solve the system with inexact \
coefficients or the system obtained by direct rationalization of \
inexact numbers present in the system. Since many of the methods used \
by Reduce require exact input, providing Reduce with an exact version \
of the system may help.



Given that rmin depends on p and q is it correct to include it in L[rmin_, p_, q_?NumericQ] :=
or should it just be L[p_, q_?NumericQ] := ?


My new code is:




Clear["Global`*"]


k[r_] = 1 + b^2/r^2

h[r_] = 1 - 1/r^4

gtt[r_] = -h[r]*k[r]^(-1)*(r/R)^2*(1 + b^2*r^2)*rh^2

gyt[r_] = -2*b^2*r^2*(k[r]^(-1)*(r/R)^2)*h[r]*rh^2

gyy[r_] = k[r]^(-1)*(r/R)^2*(1 - b^2*r^2*h[r])*rh^2

gzz[r_] = (r/R)^2*rh^2

grr[r_] = (R/r)^2*h[r]^(-1)
gxx[r_] = (r/R)^2*rh^2

s = (c^2 - 1)^(1/2)

ggtt[r_] = FullSimplify[c^2*gtt[r] + s^2*gzz[r]]
ggyt[r_] = FullSimplify[c*gyt[r]]
ggzt[r_] = FullSimplify[-2*c*s*(gtt[r] + gzz[r])]
ggyz[r_] = FullSimplify[-s*gyt[r]]
ggzz[r_] = FullSimplify[c^2*gzz[r] + s^2*gtt[r]]
ggyy[r_] = FullSimplify[K[r]^(-1)*(r/R)^2*(b*R)^2*(1 - h[r])]
ggrr[r_] = FullSimplify[(R/r)^2*h[r]^(-1)]
ggxx[r_] = FullSimplify[gxx[r]]


(* Note dil[r]=e^(-2\[CapitalPhi])*)

dil[r_] = k[r]

c = 1.0
b = 1.0
rh = 1
R = 1
\[CapitalLambda] = 100000


xsigma[r_] = p/q*(-ggtt[r]*ggzz[r] + ggzt[r]^2)/(-ggtt[r]*ggxx[r])


rsigma[r_] =
1/q/(-ggtt[r]*ggrr[r])^(0.5)*(-ggtt[r]*ggzz[r] +
ggzt[r]^2)^(0.5)*((-ggtt[r]*ggzz[r] + ggzt[r]^2) (dil[r] +
p^2/(ggtt[r]*ggxx[r])) - q^2)^(0.5)



rmin[p_, q_?NumericQ] :=
NSolve[(-ggtt[r]*ggzz[r] + ggzt[r]^2) (dil[r] +
p^2/(ggtt[r]*ggxx[r])) - q^2 == 0, r]


Theta[rmin_, p_, q_?NumericQ] :=
NIntegrate[xsigma[r]/rsigma[r], {r, rmin, \[CapitalLambda]}]/
NIntegrate[1/rsigma[r], {r, rmin, \[CapitalLambda]}]


L[rmin_, p_, q_?NumericQ] :=
2*((NIntegrate[
1/rsigma[r], {r, rmin, \[CapitalLambda]}])^2 + (NIntegrate[
xsigma[r]/rsigma[r], {r, rmin, \[CapitalLambda]}])^2)^(0.5)

S[rmin_, p_, q_?NumericQ] :=
dil[r]/q*(ggzt[r]^2 - ggtt[r]*ggzz[r])/rsigma[r]


qOfpTheta[p_, Theta_] := NSolve[Theta == Theta[p, q], q][[1, 2]]

Theta = 0.8

ParametricPlot[{L[p, qOfpTheta[p, Theta]],
S[p, qOfpTheta[p, Theta]]}, {p, 0, 10},
AspectRatio -> 1/GoldenRatio, Frame -> True, Axes -> False]

 

[Dale]

Thanks for your reply. I have made the changes you suggested, however I am getting the following error messages when I run my new code:

Inverse functions are being used by Solve, so some solutions may not \
be found; use Reduce for complete solution information.

NSolve is really only meant for polynomials, so I guess it is choking on your expression. You may need to use FindRoot, however with FindRoot be forewarned that if there are multiple solutions it will converge to the one closest to your initial guess.

Given that rmin depends on p and q is it correct to include it in L[rmin_, p_, q_?NumericQ] :=
or should it just be L[p_, q_?NumericQ] := ?

My personal preference would be the second style, otherwise you risk putting in a bad expression for rmin.

ggyy[r_] = FullSimplify[K[r]^(-1)*(r/R)^2*(b*R)^2*(1 - h[r])]

I think you mean lower case K

c = 1.0
b = 1.0
rh = 1
R = 1
\[CapitalLambda] = 100000


xsigma[r_] = p/q*(-ggtt[r]*ggzz[r] + ggzt[r]^2)/(-ggtt[r]*ggxx[r])

Did you know that if you replace c and b with exact numbers (c=1 and b=1) you get that xsigma FullSimplifies to p/q. Is that correct?

rsigma[r_] =
1/q/(-ggtt[r]*ggrr[r])^(0.5)*(-ggtt[r]*ggzz[r] +
ggzt[r]^2)^(0.5)*((-ggtt[r]*ggzz[r] + ggzt[r]^2) (dil[r] +
p^2/(ggtt[r]*ggxx[r])) - q^2)^(0.5)

rsigma also FullSimplifies considerably if you replace the 0.5 exponents with 1/2 and if you can assume 0<r

rmin[p_, q_?NumericQ] :=
NSolve[(-ggtt[r]*ggzz[r] + ggzt[r]^2) (dil[r] +
p^2/(ggtt[r]*ggxx[r])) - q^2 == 0, r]

With the above simplifications you can solve this equation for rmin reduces to a very simple 6th order polynomial in r which can be solved analytically. Since it is sixth order there are, of course, 6 roots and you will have to pick the right one, but that should still be preferable to numerically solving it. My guess is you want the second root.

qOfpTheta[p_, Theta_] := NSolve[Theta == Theta[p, q], q][[1, 2]]

Theta = 0.8

The highlighted theta needs to have a different name from the function Theta. Try lower-case theta.

 

[hufflepup]

I have implemented those changes however I am getting error messages saying that the intergrand has evaluated to non numeric values in all parts of the sampling reason, I think this is because it is not inputting the values for q correctly but I can't figure out why.

The expression does simplify for b=c=1, however I am also interested in these functions at other values of b and c, Mathematica can still find a solution in this case, but it is only valid over a range of values for p and q (which is what I expect for the behaviour of these functions).

What does the [1,2] in

qOfpTheta[p_, theta_] :=
FindRoot[theta == Theta[p, q], {q, 10}][[1, 2]]

mean?



Clear["Global`*"]


c = 2.0
b = 1.0
rh = 1
R = 1
\[CapitalLambda] = 100000


k[r_] = 1 + b^2/r^2

h[r_] = 1 - 1/r^4

gtt[r_] = -h[r]*k[r]^(-1)*(r/R)^2*(1 + b^2*r^2)*rh^2

gyt[r_] = -2*b^2*r^2*(k[r]^(-1)*(r/R)^2)*h[r]*rh^2

gyy[r_] = k[r]^(-1)*(r/R)^2*(1 - b^2*r^2*h[r])*rh^2

gzz[r_] = (r/R)^2*rh^2

grr[r_] = (R/r)^2*h[r]^(-1)
gxx[r_] = (r/R)^2*rh^2

s = (c^2 - 1)^(1/2)

ggtt[r_] = FullSimplify[c^2*gtt[r] + s^2*gzz[r]]
ggyt[r_] = FullSimplify[c*gyt[r]]
ggzt[r_] = FullSimplify[-2*c*s*(gtt[r] + gzz[r])]
ggyz[r_] = FullSimplify[-s*gyt[r]]
ggzz[r_] = FullSimplify[c^2*gzz[r] + s^2*gtt[r]]
ggyy[r_] = FullSimplify[k[r]^(-1)*(r/R)^2*(b*R)^2*(1 - h[r])]
ggrr[r_] = FullSimplify[(R/r)^2*h[r]^(-1)]
ggxx[r_] = FullSimplify[gxx[r]]


(* Note dil[r]=e^(-2\[CapitalPhi])*)

dil[r_] = k[r]



xsigma[r_] = p/q*(-ggtt[r]*ggzz[r] + ggzt[r]^2)/(-ggtt[r]*ggxx[r])


rsigma[r_] =
1/q/(-ggtt[r]*ggrr[r])^(1/2)*(-ggtt[r]*ggzz[r] +
ggzt[r]^2)^(0.5)*((-ggtt[r]*ggzz[r] + ggzt[r]^2) (dil[r] +
p^2/(ggtt[r]*ggxx[r])) - q^2)^(1/2)



rmin[p_, q_] := \[Sqrt]Root[
c^2 (3 p^2 R^4 - 3 c^2 p^2 R^4 + 3 c^2 rh^4 - 3 c^4 rh^4) +
b^2 (p^2 R^4 + 12 c^2 p^2 R^4 - 12 c^4 p^2 R^4 + c^2 q^2 R^4 -
2 c^2 rh^4 + 21 c^4 rh^4 - 18 c^6 rh^4) #1 + (p^2 R^4 +
b^4 p^2 R^4 + 12 b^4 c^2 p^2 R^4 - 12 b^4 c^4 p^2 R^4 -
b^4 q^2 R^4 + c^2 q^2 R^4 + 2 b^4 c^2 q^2 R^4 - b^4 rh^4 -
2 c^2 rh^4 - 9 b^4 c^2 rh^4 + 3 c^4 rh^4 + 48 b^4 c^4 rh^4 -
36 b^4 c^6 rh^4) #1^2 +
b^2 (-6 c^2 p^2 R^4 + 6 c^4 p^2 R^4 - 2 q^2 R^4 +
2 c^2 q^2 R^4 - 2 rh^4 - b^4 rh^4 - 10 c^2 rh^4 -
10 b^4 c^2 rh^4 + 3 c^4 rh^4 + 36 b^4 c^4 rh^4 + 9 c^6 rh^4 -
24 b^4 c^6 rh^4) #1^3 + (-p^2 R^4 - b^4 p^2 R^4 -
12 b^4 c^2 p^2 R^4 + 12 b^4 c^4 p^2 R^4 - q^2 R^4 -
b^4 c^2 q^2 R^4 - rh^4 - b^4 rh^4 - c^2 rh^4 -
8 b^4 c^2 rh^4 - 30 b^4 c^4 rh^4 + 36 b^4 c^6 rh^4) #1^4 +
b^2 (-p^2 R^4 - c^2 q^2 R^4 + rh^4 + b^4 rh^4 + 2 c^2 rh^4 +
9 b^4 c^2 rh^4 - 6 c^4 rh^4 - 48 b^4 c^4 rh^4 +
36 b^4 c^6 rh^4) #1^5 + (3 b^4 c^2 p^2 R^4 -
3 b^4 c^4 p^2 R^4 + rh^4 + 2 b^4 rh^4 + 10 b^4 c^2 rh^4 -
3 b^4 c^4 rh^4 - 9 b^4 c^6 rh^4) #1^6 +
b^2 (1 + c^2 - 2 b^4 c^2 + 21 b^4 c^4 - 18 b^4 c^6) rh^4 #1^7 +
b^4 c^2 (-2 + 3 c^2) rh^4 #1^8 +
b^6 c^4 (-3 + 3 c^2) rh^4 #1^9 &, 7]

Theta[p_, q_?NumericQ] :=
NIntegrate[xsigma[r]/rsigma[r], {r, rmin[p, q], \[CapitalLambda]}]/
NIntegrate[1/rsigma[r], {r, rmin[p, q], \[CapitalLambda]}]


L[p_, q_?NumericQ] :=
2*((NIntegrate[
1/rsigma[r], {r,
rmin[p, q], \[CapitalLambda]}])^2 + (NIntegrate[
xsigma[r]/rsigma[r], {r,
rmin[p, q], \[CapitalLambda]}])^2)^(0.5)

S[p_, q_?NumericQ] :=
NIntegrate[
dil[r]/q*(ggzt[r]^2 - ggtt[r]*ggzz[r])/rsigma[r], {r,
rmin[p, q], \[CapitalLambda]}]


qOfpTheta[p_, theta_] :=
FindRoot[theta == Theta[p, q], {q, 10}][[1, 2]]

theta = 1000

ParametricPlot[{L[p, qOfpTheta[p, theta]],
S[p, qOfpTheta[p, theta]]}, {p, 1, 2},
AspectRatio -> 1/GoldenRatio, Frame -> True, Axes -> False]