3
$\begingroup$

The code below is a simplified example of a problem I am working on. I first define a relation V that depends on x[t] and y[t], and I also define its derivative Vx. Next, I solve V==0 for the variable x[t] in terms of y[t]:

V = x[t]^3 - (x[t] + y[t])^0.5;
Vx = D[V, x[t]];
x[t] = x[t] /. Solve[V == 0, x[t]];
Vx /. y[t] -> 20

The last expression shows that the vector x[t] has 6 elements. However, I am only interested in the value of x[t] that corresponds to the smallest real value of Vx. The following commands choose this value and call it x1[t]:

Vx0 = (Vx /. x[t] -> #) &;
Vx1 = # &;
Vx2 = If[Abs[Im[Vx1[#]]] < 10^(-10), Re[Vx1[#]], 10^8] &;
Vx3 = Vx2[#] & /@ Vx;
x1[t] = Pick[x[t], Vx3, Min[Vx3]];

The last line gives me an error that the Pick[] expressions have incompatible shapes. However, Length[Vx3] = Length[x[t] = 6. Moreover, I get a numerical solution upon specifying a value of y[t]. For instance;

x1[t] /. y[t] -> 20

yields {-1.62445138605743}. I can also plot x1[t] as a function of y[t]:

Plot[{x1[t]}, {y[t], 0, 50}, PlotRange -> {-5, 5}, AspectRatio -> .5]

enter image description here

yields a nice graph. I am having two difficulties with this type of model. First, I cannot always get the Plot command to work properly in my more complex model. Second, I cannot substitute x1[t] into a differential equation to solve it numerically:

soly = NDSolve[{Derivative[1][y][t] == 
     0.05*y[t]*(1 - y[t]/50) - x1[t], y[0] == 70}, {y}, {t, 0, 20}];

does not work. Instead it tells me that I have incompatible shapes. I really thought this would work given that the numerical plot worked. Any ideas?

$\endgroup$
1
  • $\begingroup$ Welcome to Mathematica.SE, Rick! I suggest the following: 1) Take the tour and check the faqs. 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – Chris K
    May 7 '19 at 11:21
4
$\begingroup$

Root objects take a pure function and an index, and the indices are ordered so that real values come before complex values, and amongst the real values they are sorted in increasing order. Now, let's look at your x[t] (after rationalizing the exponent):

V = x^3 - (x + y)^(1/2);
x /. Solve[V == 0, x]

{Root[-y - #1 + #1^6 &, 1], Root[-y - #1 + #1^6 &, 2], Root[-y - #1 + #1^6 &, 3], Root[-y - #1 + #1^6 &, 4], Root[-y - #1 + #1^6 &, 5], Root[-y - #1 + #1^6 &, 6]}

If there is a minimum real root, it will be the first one above, so you can define:

x[t_] := Root[-y[t] - # + #^6&, 1]

Then, a plot of x[t] gives:

Plot[x[t], {y[t], 0, 50}, PlotRange->{-5,5}]

enter image description here

in agreement with your plot. Your ODE can use the above x[t]:

sol = NDSolveValue[{y'[t] == .05 y[t] (1-y[t]/50) - x[t], y[0]==70}, y, {t, 0, 20}];

Visualization:

Plot[sol[t], {t, 0, 20}]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks! This did help. $\endgroup$
    – Rick
    May 7 '19 at 14:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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