Parametric Plots in Mathematica

In summary: Theta[rmin_, p_, q_?NumericQ] := NIntegrate[xsigma[r]/rsigma[r], {r, rmin, \[CapitalLambda]}]/ NIntegrate[1/rsigma[r], {r, rmin, \[CapitalLambda]}]... can be substituted into rmin (for a given p and q) this has the potential to speed things up. It is still a bit of a mess, but at least it is a mess that can be solved with Root objects.
  • #1
hufflepup
6
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
  • #2
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]
 
  • #3
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]
 
  • #4
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.
 
  • #5
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]
 

1. What is a parametric plot in Mathematica?

A parametric plot in Mathematica is a type of plot that uses parametric equations to represent a curve or surface. Parametric equations define the coordinates of points on a curve or surface in terms of one or more parameters. These plots are useful for visualizing complex shapes and functions that cannot be easily represented with traditional Cartesian coordinates.

2. How do I create a parametric plot in Mathematica?

To create a parametric plot in Mathematica, you first need to define the parametric equations for the curve or surface you want to plot. Then, use the ParametricPlot function and input your equations as the first argument. You can also customize the appearance of the plot by specifying options such as colors, axes labels, and plot ranges.

3. Can I plot multiple curves on the same parametric plot?

Yes, you can plot multiple curves on the same parametric plot by using the ParametricPlot function with a list of parametric equations as the first argument. Each equation will be plotted with a different color by default, but you can also specify custom colors for each curve.

4. How can I incorporate parameters in my parametric plot?

To incorporate parameters in your parametric plot, simply include them in your parametric equations. You can then manipulate these parameters to see how they affect the shape of the curve or surface. This is useful for exploring the behavior of functions and creating interactive visualizations.

5. Are there any other useful functions for working with parametric plots in Mathematica?

Yes, Mathematica has several other functions that can be useful for working with parametric plots. These include ParametricRegion for creating regions defined by parametric equations, ParametricPlot3D for plotting surfaces in 3D, and ParametricNDSolve for solving differential equations with parameters. You can also use the Manipulate function to create interactive parametric plots that allow you to change the values of parameters in real-time.

Similar threads

Every integer n>1 is the product of a square-free integer?
    • Calculus and Beyond Homework Help
Replies
5
Views
2K
MATLAB How can I plot a Hohmann Transfer Orbit in MATLAB using ode45?
    • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
Action variables for the Kepler problem
    • Introductory Physics Homework Help
Replies
2
Views
1K
Ideal gas in a cylindrical container
    • Advanced Physics Homework Help
Replies
4
Views
1K
Different sunset times due to elevation ##h## at a point on the Earth
    • Introductory Physics Homework Help
Replies
13
Views
489
I Atmospheric Density vs Altitude in an O'Neil Cylinder
    • Classical Physics
Replies
1
Views
587
Plotting points in three-dimensional space
    • Precalculus Mathematics Homework Help
Replies
8
Views
884
Simple Induction Induction proof of Polynomial Division Theorem
    • Math Proof Training and Practice
Replies
10
Views
1K
Mathematica Verifying PDE solutions using Mathematica
    • MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
2K
Mathematica Mathematica - using NIntegrate within NDSolve
    • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
3K

Hot Threads

Recent Insights

Back
Top