Parametric Plots in Mathematica

Click For Summary

Discussion Overview

The discussion revolves around the implementation of parametric plots in Mathematica, specifically focusing on plotting the function S against L for given parameters p and q, while keeping Theta fixed. Participants explore numerical solutions, code corrections, and potential errors in the provided Mathematica code.

Discussion Character

  • Technical explanation
  • Debate/contested
  • Mathematical reasoning

Main Points Raised

  • One participant seeks to plot S against L for fixed Theta using Mathematica, mentioning the complexity of rmin's dependence on p and q.
  • Another participant suggests renaming the symbol K to avoid conflicts with built-in symbols and points out undefined variables gtx and gtz in the code.
  • Concerns are raised about the use of NSolve versus Reduce for solving equations, with a recommendation to use FindRoot for numerical solutions.
  • There are discussions about the implications of including rmin in the function definitions for L and S, with differing preferences expressed.
  • Participants note that certain expressions simplify significantly under specific conditions, which could lead to easier analytical solutions.
  • Errors related to inverse functions and the need for exact coefficients are highlighted, indicating potential issues with the numerical methods employed.
  • One participant points out that the naming of the variable Theta conflicts with the function Theta, suggesting a need for distinct naming conventions.

Areas of Agreement / Disagreement

Participants express differing opinions on the best approach to solving the equations and the structure of the code. There is no consensus on the optimal method for implementing the parametric plot or resolving the encountered errors.

Contextual Notes

Limitations include unresolved mathematical steps, the dependence on specific definitions, and the potential for multiple solutions in numerical methods. The discussion reflects ongoing refinements and corrections to the code without reaching definitive conclusions.

hufflepup
Messages
6
Reaction score
0
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?
 
Physics news on Phys.org
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]
 
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]
 
hufflepup said:
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.
hufflepup said:
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.

hufflepup said:
ggyy[r_] = FullSimplify[K[r]^(-1)*(r/R)^2*(b*R)^2*(1 - h[r])]
I think you mean lower case K
hufflepup said:
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?

hufflepup said:
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
hufflepup said:
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.


hufflepup said:
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.
 
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]
 

Similar threads

Mathematica Coloring and legending issue
  • · Replies 1 ·
Replies
1
Views
3K
Mathematica Solving equation and plotting issue with Mathematica
  • · Replies 2 ·
Replies
2
Views
3K
Mathematica Solve equation perturbatively
  • · Replies 2 ·
Replies
2
Views
2K
Every integer n>1 is the product of a square-free integer?
  • · Replies 5 ·
Replies
5
Views
4K
Mathematica Graph Plotting Problem in Mathematica
  • · Replies 2 ·
Replies
2
Views
2K
Mathematica Rotation of 3D Plot using Euler angles
  • · Replies 2 ·
Replies
2
Views
3K
Mathematica Plotting Mathematica f[r_][x_]: Solving Derivatives & Plotting Lines
  • · Replies 1 ·
Replies
1
Views
2K
Action variables for the Kepler problem
  • · Replies 2 ·
Replies
2
Views
3K
Ideal gas in a cylindrical container
  • · Replies 4 ·
Replies
4
Views
2K
Mathematica Verifying PDE solutions using Mathematica
  • · Replies 5 ·
Replies
5
Views
3K