3
$\begingroup$

I want to create a simple code to solve an algebraic equation named eq[x], and plot x as a function of three parameters y, z, and k (not with Plot3D). To do this, I fix y and z each time and plot x as a function of k: y will take the values ​​2,5,20, while z will take 0.1,1,10,100,200, with k varying between 0 and 50.

I'd like to achieve all this in a single subroutine without having to systematically use the Clear["Global"] instruction to reset the values ​​of y z,k,root[k].

What's the best numerical scheme for this problem?

Here is the code that I usually use to solve and plot this equation:

eq[x_?NumericQ, k_?NumericQ] = 1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y);
Solution for y=2 and z = 0.1
 Clear["Global`*"] (*============== Solution for y=2 ==============*) y = 2; (*== z=0.1 ==*) z = 0.1; root[k_] :=SolveValues[1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, x] (*Tables of the 3 roots with k between 0 and 50 *) x1 = Table[Evaluate[{k, root[k][[1]]}], {k, 0.1, 50, 0.1}]; x2 = Table[Evaluate[{k, root[k][[2]]}], {k, 0.1, 50, 0.1}]; x3 = Table[Evaluate[{k, root[k][[3]]}], {k, 0.1, 50, 0.1}]; (*Plots of the 3 roots as a function of k*) ListLinePlot[{x1, x2, x3}, PlotRange -> All, Frame -> True, FrameLabel -> {k, x}, PlotLegends -> {"x1,", "x2", "x3"}]
Solution for y=2 and z = 1
Clear["Global`*"] (*============== Solution for y=2 ==============*) y = 2; (*== z=1 ==*) z = 1; root[k_] :=SolveValues[1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, x] (*Tables of the 3 roots with k between 0 and 50 *) x1 = Table[Evaluate[{k, root[k][[1]]}], {k, 0.1, 50, 0.1}]; x2 = Table[Evaluate[{k, root[k][[2]]}], {k, 0.1, 50, 0.1}]; x3 = Table[Evaluate[{k, root[k][[3]]}], {k, 0.1, 50, 0.1}]; (*Plots of the 3 roots as a function of k*) ListLinePlot[{x1, x2, x3}, PlotRange -> All, Frame -> True, FrameLabel -> {k, x}, PlotLegends -> {"x1,", "x2", "x3"}]
Solution for y=2 and z = 10
Clear["Global`*"] (*============== Solution for y=2 ==============*) y = 2; (*== z=10 ==*) z = 10; root[k_] :=SolveValues[1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, x] (*Tables of the 3 roots with k between 0 and 50 *) x1 = Table[Evaluate[{k, root[k][[1]]}], {k, 0.1, 50, 0.1}]; x2 = Table[Evaluate[{k, root[k][[2]]}], {k, 0.1, 50, 0.1}]; x3 = Table[Evaluate[{k, root[k][[3]]}], {k, 0.1, 50, 0.1}]; (*Plots of the 3 roots as a function of k*) ListLinePlot[{x1, x2, x3}, PlotRange -> All, Frame -> True, FrameLabel -> {k, x}, PlotLegends -> {"x1,", "x2", "x3"}]
Solution for y=5 and z = 0.1
 Clear["Global`*"] (*============== Solution for y=5 ==============*) y = 5; (*== z=0.1 ==*) z = 0.1; root[k_] :=SolveValues[1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, x] (*Tables of the 3 roots with k between 0 and 50 *) x1 = Table[Evaluate[{k, root[k][[1]]}], {k, 0.1, 50, 0.1}]; x2 = Table[Evaluate[{k, root[k][[2]]}], {k, 0.1, 50, 0.1}]; x3 = Table[Evaluate[{k, root[k][[3]]}], {k, 0.1, 50, 0.1}]; (*Plots of the 3 roots as a function of k*) ListLinePlot[{x1, x2, x3}, PlotRange -> All, Frame -> True, FrameLabel -> {k, x}, PlotLegends -> {"x1,", "x2", "x3"}]
Solution for y=5 and z = 1
 Clear["Global`*"] (*============== Solution for y=5 ==============*) y = 5; (*== z=1 ==*) z =1; root[k_] :=SolveValues[1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, x] (*Tables of the 3 roots with k between 0 and 50 *) x1 = Table[Evaluate[{k, root[k][[1]]}], {k, 0.1, 50, 0.1}]; x2 = Table[Evaluate[{k, root[k][[2]]}], {k, 0.1, 50, 0.1}]; x3 = Table[Evaluate[{k, root[k][[3]]}], {k, 0.1, 50, 0.1}]; (*Plots of the 3 roots as a function of k*) ListLinePlot[{x1, x2, x3}, PlotRange -> All, Frame -> True, FrameLabel -> {k, x}, PlotLegends -> {"x1,", "x2", "x3"}]

.... and so on.

Improve this question
$\endgroup$

    4 Answers 4

    Reset to default
    6
    $\begingroup$

    Your equation is a fifth-order polynomial and Mathematica can represent its solutions as Root objects:

    SolveValues[1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, x] (* Root[..., 1], Root[..., 2], Root[..., 3], Root[..., 4], Root[..., 5] *)

    Note that the Abel–Ruffini theorem proves that there is no easier way to write these roots.

    We can define a function that evaluates the $n^{\text{th}}$ such root: taking the above roots and simplifying their argument,

    R[k_, y_, z_, n_] := Root[-k^2 y^3 # - 5 (# - 1)^2 (z + BesselI[3, z]) (# y^2 BesselI[2, z] + 3 (k^2 - #^2) (# - 1) BesselI[3, z]) &, n]

    we can plot these directly, without making a table of values and without having to avoid the $k=0$ point:

    Plot[Evaluate@Table[R[k, 2, 0.1, n], {n, 5}], {k, 0, 50}]

    enter image description here

    Roots $n=4$ and $n=5$ are complex and don't show up in this plot. Use ReImPlot instead of Plot to show these complex roots too.

    This can be plotted in any old way and nothing needs to be cleared between plots. For example, making a grid of the plots you want for the given lists of $(y,z)$ pairs:

    GraphicsGrid[Table[Plot[Evaluate@Table[R[k, y, z, n], {n, 5}], {k, 0, 50}], {y, {2, 5, 20}}, {z, {0.1, 1, 10, 100, 200}}]]

    enter image description here

    Improve this answer
    $\endgroup$
    1
    • $\begingroup$Excellent @Roman thank you so much.$\endgroup$
      – Gallagher
      Commented14 hours ago
    3
    $\begingroup$

    We simply pack your code into a function: plotSol, like:

    plotSol[y_, z_] := Module[{x1, x2, x3, root}, root[k_] := SolveValues[ 1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/ 3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, x]; x1 = Table[Evaluate[{k, root[k][[1]]}], {k, 0.1, 50, 0.1}]; x2 = Table[Evaluate[{k, root[k][[2]]}], {k, 0.1, 50, 0.1}]; x3 = Table[Evaluate[{k, root[k][[3]]}], {k, 0.1, 50, 0.1}]; (*Plots of the 3 roots as a function of k*) ListLinePlot[{x1, x2, x3}, PlotRange -> All, Frame -> True, FrameLabel -> {k, x}, PlotLegends -> {"x1,", "x2", "x3"}] ]

    With this:

    plotSol[2, 0.1]

    enter image description here

    plotSol[2, 1]

    enter image description here

    plotSol[2, 10]

    enter image description here

    e.t.c.

    Improve this answer
    $\endgroup$
    1
    • $\begingroup$Many thanks @Daniel$\endgroup$
      – Gallagher
      Commented14 hours ago
    3
    $\begingroup$
    $Version (* "14.2.1 for Mac OS X ARM (64-bit) (March 16, 2025)" *) Clear["Global`*"]

    Memoize roots to speed up Manipulate

    root[k_, y_, z_] := root[k, y, z] = SolveValues[{1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/ 3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0, 0 < k < 50}, x, Reals] Table[root[k, y, z], {y, {2, 5}}, {z, {1/10, 1, 10}}]; Manipulate[ Plot[Evaluate@root[k, y, z], {k, 0, 50}, Frame -> True, FrameLabel -> (Style[#, 14] & /@ {k, root}), PlotLegends -> (Subscript["x", #] & /@ Range[3]), PlotPoints -> 50, MaxRecursion -> 5, WorkingPrecision -> 15], Row[{ Control[{{y, 2}, {2, 5}}], Spacer[40], Control[{{z, 1/10}, {1/10, 1, 10}}]}], SynchronousUpdating -> False, TrackedSymbols :> {y, z}]

    enter image description here

    EDIT: Without Manipulate

    Legended[ Table[ Plot[Evaluate@root[k, y, z], {k, 0, 50}, Frame -> True, FrameLabel -> (Style[#, 14] & /@ {k, root}), PlotLabel -> StringForm["y = ``; z = ``", y, z], PlotPoints -> 50, MaxRecursion -> 5, WorkingPrecision -> 15], {y, {2, 5}}, {z, {1/10, 1, 10}}] // Grid, LineLegend[ColorData[97] /@ Range[3], (Subscript["x", #] & /@ Range[3])]]

    enter image description here

    Improve this answer
    $\endgroup$
    5
    • $\begingroup$Thank you dear @Bob. So without Manipulate what can we do?$\endgroup$
      – Gallagher
      Commentedyesterday
    • $\begingroup$Thank you so much @Bob. Usually I solve more complicated equations. In this case how to proceed with FindRoot instead of SolveValues?$\endgroup$
      – Gallagher
      Commentedyesterday
    • 2
      $\begingroup$Don’t ask new questions in comments. Post a new question with a specific example of the problem you are experiencing.$\endgroup$
      – Bob Hanlon
      Commentedyesterday
    • $\begingroup$Bob how to proceed with FindRoot with the same equation in this post.$\endgroup$
      – Gallagher
      Commentedyesterday
    • 2
      $\begingroup$Post a new question, show what you have tried, and explain with what you are having a problem.$\endgroup$
      – Bob Hanlon
      Commentedyesterday
    1
    $\begingroup$

    version ContourPlot

    yi = {2, 5, 20}; zi = {0.1 , 1, 10, 100, 200 }; sol[y_ ] := Show@Map[ContourPlot[ Evaluate[ 1 - k^2/x^2 - y^2/(BesselI[3, z] x^2 (1 - 1/x)) (BesselI[2, z]/ 3 + (k^2)/((BesselI[3, z] + z) 15 (1 - 1/x)^2 x^2) y) == 0 /. z -> # ], {k, 0, 50}, {x, -300, 300}, ContourStyle -> RandomColor[], PlotPoints -> 100, PlotLegends -> {"y=" <> ToString[y] <> " z=" <> ToString[#] } ] &, zi] Column[Map[sol, {2, 5, 20}] ]

    enter image description here

    Improve this answer
    $\endgroup$
    2
    • $\begingroup$Thank you so much @Neumann$\endgroup$
      – Gallagher
      Commented14 hours ago
    • $\begingroup$You are welcome$\endgroup$
      – Ulrich Neumann
      Commented14 hours ago

    Start asking to get answers

    Find the answer to your question by asking.

    Ask question

    Explore related questions

    See similar questions with these tags.