Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Parametric Plots in Mathematica

  1. Dec 18, 2009 #1
    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?
     
  2. jcsd
  3. Dec 19, 2009 #2

    Dale

    Staff: Mentor

    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]
     
  4. Dec 19, 2009 #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]
     
  5. Dec 19, 2009 #4

    Dale

    Staff: Mentor

    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.
    My personal preference would be the second style, otherwise you risk putting in a bad expression for rmin.

    I think you mean lower case K
    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 also FullSimplifies considerably if you replace the 0.5 exponents with 1/2 and if you can assume 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.


    The highlighted theta needs to have a different name from the function Theta. Try lower-case theta.
     
  6. Dec 20, 2009 #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]
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook