How to plot a function with two parameters determined by another parameter?

Question

I have one small question. My code is like follows.

a = 5; 
mm = Solve[a x + y == 7 && x - y == 1, {x, y}]; 
t1 = x /. mm[[1]]; 
t2 = y /. mm[[1]]; 

myfunction[{x_, y_}] := (y + x)/(x + r)
Plot[myfunction /@ {{t1, t2}}, {r, 0, 3}, Evaluated -> True]

I want to change plot of myfunction for different a in one plot. How could I do? Thanks.

Edit - The actual code is more complex than the stripped down example above and fails to produce the figure showing the function ppw[r] for different values of bb. I have tried the replacement methods suggested by Peltio, but they do not work

Here is a simplified version of the code using only linear ODEs:

xxin[aa_?NumericQ, bb_?NumericQ] :=NDSolve[{\[Psi]''[r] == -bb \[Psi][r], \[Psi]'[0.1] ==0, \[Psi][1] == aa}, \[Psi], {r, 0.1, 1}, PrecisionGoal -> 20]
xxshell[aa_?NumericQ, cc_?NumericQ] :=NDSolve[{\[Psi]''[r] == \[Psi][r] + 1, \[Psi][1] ==aa, \[Psi][1.5] == cc}, \[Psi], {r, 1, 1.5}, PrecisionGoal -> 20]
xxout[cc_?NumericQ] :=NDSolve[{\[Psi]''[r] == \[Psi][r], \[Psi][1.5] ==cc, \[Psi]'[3] == 0}, \[Psi], {r, 1.5, 3}, PrecisionGoal -> 20]
xxinaTry[aa_?NumericQ,bb_?NumericQ] := (Evaluate[\[Psi]'[1] /. xxin[aa, bb]])[[1]]
xxoutaTry[aa_?NumericQ,cc_?NumericQ] := (Evaluate[\[Psi]'[1] /. xxshell[aa, cc]])[[1]]
xxinbTry[aa_?NumericQ,cc_?NumericQ] := (Evaluate[\[Psi]'[1.5] /. xxshell[aa, cc]])[[1]]
xxoutbTry[cc_?NumericQ] := (Evaluate[\[Psi]'[1.5] /. xxout[cc]])[[1]]

 bb = 0.8; 
 sols = FindRoot[{ xxinaTry[aa, bb] == bb xxoutaTry[aa, cc], 
 xxoutbTry[cc] == xxinbTry[aa, cc]}, {{aa, 0}, {cc, 0}}]; 
\[Psi]a =aa /. sols;
\[Psi]b = cc /. sols;
pww[r_] = Evaluate[Piecewise[{{Evaluate[ Exp[\[Psi][r]] /. xxin[\[Psi]a, bb]], 
 0.1 <= r <= 1}, {Evaluate[
  Exp[\[Psi][r]] /. xxshell[\[Psi]a, \[Psi]b]], 
 1 < r < 1.5}, {Evaluate[Exp[\[Psi][r]] /. xxout[\[Psi]b]], 
 r > 1.5}}, Indeterminate]]; 
 Show[Plot[pww[r], {r, 0, 3}], PlotRange -> All]

Edit 2: By using the above simplified code, I successfully created a plot in the end. (sorry, I do not know how to display a figure here)

But when I modified the pww[r] function into the following:

 pww[r_][bb_?NumericQ] = Evaluate[Piecewise[{{Evaluate[ Exp[\[Psi][r]] /. xxin[\[Psi]a, bb]], 
 0.1 <= r <= 1}, {Evaluate[
  Exp[\[Psi][r]] /. xxshell[\[Psi]a, \[Psi]b]], 
 1 < r < 1.5}, {Evaluate[Exp[\[Psi][r]] /. xxout[\[Psi]b]], 
 r > 1.5}}, Indeterminate]]; 
  Show[Table[Plot[pww[r][bb], {r, 0, 3}], {bb, 0.8, 0.8, 0.1}], PlotRange -> All]

Here, bb is a variable and I delete the bb=0.8; Then

   sols = FindRoot[{ xxinaTry[aa, bb] == bb xxoutaTry[aa, cc], xxoutbTry[cc] == xxinbTry[aa, cc]}, {{aa, 0}, {cc, 0}}];

does not work because it does not know bb. In this case, I tried to replace the \[Psi]a and \[Psi]b in the function pww[r_][bb_?NumericQ] by using the following code:

   Evaluate[aa /. FindRoot[{xxinaTry[aa, bb] == bb xxoutaTry[aa, cc], 
xxoutbTry[cc] == xxinbTry[aa, cc]}, {{aa, 0}, {cc, 0}}]]
   Evaluate[cc /. FindRoot[{xxinaTry[aa, bb] == bb xxoutaTry[aa, cc], 
xxoutbTry[cc] == xxinbTry[aa, cc]}, {{aa, 0}, {cc, 0}}]]

In this case, I got the following kind of errors:

During evaluation of In[1]:= FindRoot::nlnum: The function value {0.244919 bb+xxinaTry[0.,bb],-0.244919} is not a list of numbers with dimensions {2} at {aa,cc} = {0.,0.}. >>

During evaluation of In[1]:= ReplaceAll::reps: {FindRoot[{xxinaTry[aa,bb]==bb xxoutaTry[aa,cc],xxoutbTry[cc]==xxinbTry[aa,cc]},{{aa,0},{cc,0}}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

During evaluation of In[1]:= NDSolve::dsvar: 0.12251026530612244` cannot be used as a variable. >>

During evaluation of In[1]:= ReplaceAll::reps: {NDSolve[{([Psi]^[Prime][Prime])[0.12251]==-0.8 [Psi][0.12251],([Psi]^[Prime])[0.1]==0,[Psi][1]==0.},[Psi],{0.12251,0.1,1},PrecisionGoal->20]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

During evaluation of In[1]:= FindRoot::nlnum: The function value {0.2 ([Psi]^[Prime])[1.],0.} is not a list of numbers with dimensions {2} at {aa,cc} = {0.,0.}. >>

...

It should be very complex. Could you provide some help? The above code can be run in Mathematica 9.0 and generate the same error in your computer, I believe.

Answer 1

Whenever you want to change something and see its effect, then Manipulate is your friend

Manipulate[

Module[{sol, x, y, r},
  sol = First@Solve[a x + y == 7 && x - y == 1, {x, y}];
  Plot[myfunction[{x /. sol, y /. sol, r}], {r, 0, 3}, Frame -> True, 
   ImageSize -> 350, ImagePadding -> 30, GridLines -> Automatic, 
   GridLinesStyle -> LightGray, FrameLabel -> {{"r", None}, {"x", "my plot"}}]
  ],

 {{a, 5, "a"}, 1, 10, .1, Appearance -> "Labeled"},

 Initialization :>
  (
   myfunction[{x_, y_, r_}] := (y + x)/(x + r) 
   )
 ]

Add

If I want to plot one figure of my function that contains a=1,3,5 for comparison, for example, how can i do this?

You can add all the values of a you want to plot to a list, then hit a button to plot them all. Move the slider to select a new value of a, then click the add button to add it. Keep doing this until you selected all the values a. Adding legend is left as an exercise to the reader ;)

may be something like this

Manipulate[tick;

 Module[{r},
  r = makeOnePlot[#] & /@ collection;
  Show[r, PlotRange -> All]
  ],

 {{a, 5, "a"}, 1, 10, .1, Appearance -> "Labeled"},
 Grid[{
   {Dynamic@collection, SpanFromLeft},
   {Button["Add", AppendTo[collection, a]],
    Button["Reset", collection = {5}; tick = Not[tick]],
    Button["Run", tick = Not[tick]]
    }}, Frame -> All, Spacings -> {1, 1}
  ],

 {{collection, {5}}, None},
 {{tick, True}, None},

 TrackedSymbols :> {tick},

 Initialization :>
  (
   myfunction[{x_, y_, r_}] := (y + x)/(x + r);
   makeOnePlot[a_] := Module[{sol, x, y, r},
     sol = First@Solve[a x + y == 7 && x - y == 1, {x, y}];
     Plot[Evaluate@myfunction[{x /. sol, y /. sol, r}], {r, 0, 3}, Frame -> True, 
      ImageSize -> 350, ImagePadding -> 30, GridLines -> Automatic, 
      GridLinesStyle -> LightGray, 
      FrameLabel -> {{"r", None}, {"x", Row[{"my plot for a=", a}]}}]
     ]

   )
 ]

And if you want to do this without Manipulate, then

myfunction[{x_, y_, r_}] := (y + x)/(x + r);
makeOnePlot[a_] := Module[{sol, x, y, r},
   sol = First@Solve[a x + y == 7 && x - y == 1, {x, y}];
   Plot[Evaluate@myfunction[{x /. sol, y /. sol, r}], {r, 0, 3}, Frame -> True, 
    ImageSize -> 400, ImagePadding -> 20, 
    GridLines -> Automatic, GridLinesStyle -> LightGray, 
    FrameLabel -> {{"r", None}, {"x", None}}]
   ];
collection = {1, 3, 5};
Labeled[Show[makeOnePlot[#] & /@ collection], collection, Top]

Answer 2

I'd go for a minimalistic approach. I'd define a function that returns the solutions for a given value of the parameter a

sols[a_] := Solve[a x + y == 7 && x - y == 1, {x, y}][[1]]

and then I'd use these solutions to define a parametric function of x and y (with parameter a). EDIT: Whoops!, in my previous - too hasty - answer, I had used x and y, while it's r and f[r]. This is the correct code:

myfunction[r_][a_] := (y + x)/(x + r) /. sols[a]

That's it. Now you can either plot your solutions in the r-f[r] plane with parameter a

Show[Table[Plot[myfunction[r][a], {r, 0, 3}], {a, 0, 2, .2}]]

or, if you parameter a can vary continuously, you could even show everything on a 3D plot

Plot3D[myfunction[r][a], {r, 0, 3}, {a, 0, 2}]

EDIT: Since I have now dropped x and y, it'd be better to make those variables local. It can be done by compacting the function code in a single Block of code, e.g.:

myfunction[r_][a_] := Block[{x, y},
    (y + x)/(x + r) /. Solve[a x + y == 7 && x - y == 1, {x, y}][[1]]
]

And then you can plot as before. Or you can generate a list of functions for the values of a you need either with a Table command or by mapping myfunction on the list of values like this

myfunList = myfunction[r][#] & /@ {1, 3, 5}

    (* {7/(4 + r), 3/(2 + r), 5/(3 (4/3 + r))} *)

Plot[Evaluate@myfunList, {r,0,3}]

Evaluate might not be necessary, but I put it in there as an habit.

EDIT 2 - It is also possible to use FindRoot by incorporating it in the procedure myfunction

myfunction[r_][a_] := Block[{x, y},
    (y + x)/(x + r) /. FindRoot[{ Exp[a x - 2] == y, y^2 == x}, {{x, 1}, {y, 1}}]
    ]

Plotting can be done as above, for a finite list of selected values of a (by mapping on the list), for a range of values in a 2D plot, or for a continuous range in a 3D plot.

Answer 3

I prefer to add another answer since I believe this post should be split into two separate posts. To answer the question posed in your second EDIT, I just noticed that in the code that gives you errors, you used Set (=) instead of SetDelayed (:=) in the last two blocks of code. Also (and more importantly), you should use functions of bb if you want your solutions to depend on it. I also eliminated a lot of unnecessary Evaluates

Let's try again.

xxin[aa_?NumericQ, bb_?NumericQ] := NDSolve[{
      ψ''[r] == -bb ψ[r],
      ψ'[0.1] == 0, ψ[1] == aa}, ψ, {r, 0.1, 1}, PrecisionGoal -> 20]
xxshell[aa_?NumericQ, cc_?NumericQ] := NDSolve[{
      ψ''[r] == ψ[r] + 1,
      ψ[1] == aa, ψ[1.5] == cc}, ψ, {r, 1, 1.5}, PrecisionGoal -> 20]
xxout[cc_?NumericQ] := NDSolve[{
      ψ''[r] == ψ[r],
      ψ[1.5] == cc, ψ'[3] == 0}, ψ, {r, 1.5, 3}, PrecisionGoal -> 20]

xxinaTry[aa_?NumericQ, bb_?NumericQ] := ψ'[1] /. xxin[aa, bb][[1]]
xxoutaTry[aa_?NumericQ, cc_?NumericQ] := ψ'[1] /. xxshell[aa, cc][[1]]
xxinbTry[aa_?NumericQ, cc_?NumericQ] := ψ'[1.5] /. xxshell[aa, cc][[1]]
xxoutbTry[cc_?NumericQ] := ψ'[1.5] /. xxout[cc][[1]]

sols[bb_] := FindRoot[{
        xxinaTry[aa, bb] == bb xxoutaTry[aa, cc],
        xxoutbTry[cc] == xxinbTry[aa, cc]},
      {{aa, 0}, {cc, 0}}
      ];

ψa[bb_] := aa /. sols[bb]; ψb[bb_] := cc /. sols[bb];

pww[bb_] := Piecewise[{
        {Exp[ψ[r] /. xxin[ψa[bb], bb]], 0.1 ≤ r ≤ 1},
        {Exp[ψ[r] /. xxshell[ψa[bb], ψb[bb]]], 1 < r < 1.5},
        {Exp[ψ[r] /. xxout[ψb[bb]]], r > 1.5}
        }, Indeterminate
      ];

Table[Plot[Evaluate[pww[bb]], {r, 0.1, 3}], {bb, 0.6, 1.1, .05}]

Working?