0
$\begingroup$

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.

$\endgroup$
0

3 Answers 3

4
$\begingroup$

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) 
   )
 ]

enter image description here

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}]}}]
     ]

   )
 ]

enter image description here

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]

Mathematica graphics

$\endgroup$
2
  • $\begingroup$ Thanks for your answers. This is very interesting. 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? I am a new user for mathematica. Thanks a lot. $\endgroup$
    – user12307
    Feb 11, 2014 at 6:13
  • $\begingroup$ Thanks very much. Your answers are very beautiful, I would like to learn it carefully. $\endgroup$
    – user12307
    Feb 11, 2014 at 15:31
2
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ Thanks for these answers. I have mapped your procedure to my questions. In fact, I have still some questions. In my questions, the equations cannot use Solve like this: sols[a_] := Solve[a x + y == 7 && x - y == 1, {x, y}][[1]]. In fact, i have to find a root such as: FindRoot[{Exp[ax - 2] == y, y^2 == x}, {{x, 1}, {y, 1}}], in this case, how to write function sols[a_] to be used in this procedure? Thanks very much $\endgroup$
    – user12307
    Feb 11, 2014 at 15:25
  • $\begingroup$ Please see my second EDIT. Could it be that you wrote ax instead of a x? $\endgroup$
    – Peltio
    Feb 11, 2014 at 16:37
  • $\begingroup$ Thanks for your help. I would like to show my original code. Please see the EDIT of my questions. Thanks $\endgroup$
    – user12307
    Feb 11, 2014 at 20:08
  • $\begingroup$ Woah! Can you get solutions from those differential equations systems, in the first place? $\endgroup$
    – Peltio
    Feb 12, 2014 at 6:31
  • $\begingroup$ Yes, these equations can generate a plot in my mathematica 9.0. But when I modified them according to your procedure, it does not work. The equation: sols = FindRoot[{ xxinaTry[aa, bb] == bb xxoutaTry[aa, cc], xxoutbTry[cc] == xxinbTry[aa, cc]}, {{aa, 0}, {cc, 0}}]; can not be written as a function of bb. Thus i do not know how to do this. Of course, i can choose different bb and then plot each one, then show all the figures together. But I guess there should be one way to plot one figure containing different values of bb by using a function. If you can do this, I would be very appreciated. $\endgroup$
    – user12307
    Feb 12, 2014 at 7:07
1
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ From your silence I can infer that it did work. $\endgroup$
    – Peltio
    Feb 20, 2014 at 10:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.