3
$\begingroup$

I have a system of differential equations :

$\frac{dX}{dt} = (0.8+0.4c_{1})X-(0.007+0.006c_{2})XY$ and

$\frac{dY}{dt} = (0.048+0.004c_{3})XY-(0.4+0.2c_{4})Y$

where $0\le c_{i}\le1$.

Varying $c_{i}$ gives me an envelope of solutions for both $X(t)$ and $Y(t)$. My goal is to find the combination of $c_{i}$ which gives the upper and lower limits of my solution.

Taking $c_{i}$ as discrete points between 0 and 1 with step size 0.1, I used the following code to obtain a solution, and managed to plot it.

sol = ParametricNDSolve[{x'[t] == (0.8 + 0.4 c1) x[t] - (0.007 + 0.006 c2) x[t] y[t],y'[t] == (0.0048 + 0.0004 c3) x[t] y[t] - (0.4 + 0.2 c4) y[t], x[0] == 200, y[0] == 80}, {x, y}, {t, 0, 50}, {c1, c2, c3, c4}]
xt = Flatten[Evaluate[Table[x[c1, c2, c3, c4][t] /. sol, {c1, 0, 1, 0.1}, {c2, 0, 1,0.1}, {c3, 0, 1, 0.1}, {c4, 0, 1, 0.1}]]];
mat = Table[Table[Evaluate[xt[[a]]], {t, 0, 25, 0.01}], {a,Range[Length[xt]]}];
l = Min /@ Transpose[mat];
m = Max /@ Transpose[mat];
ind = Range[0, 25, 0.01];
lupd = Transpose[{ind, l}];
mupd = Transpose[{ind, m}];
ListLinePlot[{lupd, mupd},GridLines -> Automatic, ImageSize -> {600, 300}]

The plot for $X(t)$ came as plot The plot for $Y(t)$ came as plot Edit 3

It is not possible to calculate the equation of the upper curve of the envelope obtained by varying $c_{i}$s. Hence, I reworked my question.

Consider the following code fragment

Plot[Evaluate@({x[1, 0, 0, 1][t], x[0, 1, 1, 1][t], x[0, 0, 0, 0][t], 
     x[0.8, 0, 0, 0.2][t]} /. sol), {t, 0, 25}]

This gives 4 different trajectories for $X(t)$. enter image description here

Consider the line $t=12$. As it cuts the enevelope, we get a minimum point and a maximum point.

I would like the combination of $c_{i}$ that give this maximum and minimum points for $t=1,2,3,...,25$.

In my solution I used Max and Min to extract these maximum and minimum points and then plotted them. I could not get the combination of $c_{i}$ for each of these points.

$\endgroup$
  • $\begingroup$ Any optimization does require that you construct a merit function; that is, what exactly are you optimizing? What should be "big", and what should be "small"? $\endgroup$ – J. M.'s discontentment Jan 24 '17 at 16:53
  • $\begingroup$ The solution of the coupled equation. $X(t)$ and $Y(t)$ in this case. $\endgroup$ – Aritra Saha Jan 24 '17 at 16:59
  • 1
    $\begingroup$ How can a time-varying oscillating solution be "small"? Isn't this the lower bounds (blue) curves that you are showing on the plots? $\endgroup$ – gpap Jan 24 '17 at 17:04
  • 1
    $\begingroup$ In short, you effectively want to maximize $\|X(t)\|_\infty$ and $\|Y(t)\|_\infty$ in the "upper envelope" case? $\endgroup$ – J. M.'s discontentment Jan 24 '17 at 17:14
  • 1
    $\begingroup$ It is possible that the $c_i$ that maximize $\|X(t)\|_\infty$ can be different from the $c_i$ that maximize $\|Y(t)\|_\infty$; what do you propose to do in that case? $\endgroup$ – J. M.'s discontentment Jan 24 '17 at 17:20
3
$\begingroup$

As I pointed out in the comments, the $c_i$ do not retain the frequency of the solutions constant so the solutions cross each other and there isn't one collection of $ c_i $ that satisfies the "maximum envelope" criterion. For different integer times it's an easy application of NMaximize:

Table[NMaximize[{x[c1, c2, c3, c4][t] /. sol, 0 <= c1 <= 1, 0 <= c2 <= 1, 0 <= c3 <= 1, 0 <= c4 <= 1}, {c1, c2, c3, c4}], {t, 1, 25}]

which produces a solution (but doesn't always converge properly so maybe some playing around with the parameters of NMaximize is necessary here):

{{311.231, {c1 -> 1., c2 -> 0., c3 -> 0., 
   c4 -> 1.}}, {195.645, {c1 -> 0.851284, c2 -> 0., c3 -> 0., 
   c4 -> 1.}}, {101.746, {c1 -> 0., c2 -> 0., c3 -> 0., 
   c4 -> 1.}}, {72.4078, {c1 -> 0.554206, c2 -> 1., c3 -> 0., 
   c4 -> 1.}}, {86.7196, {c1 -> 0.516399, c2 -> 1., c3 -> 0., 
   c4 -> 1.}}, {119.434, {c1 -> 0.783024, c2 -> 1., c3 -> 0., 
   c4 -> 1.}}, {174.87, {c1 -> 1., c2 -> 1., c3 -> 0., 
   c4 -> 1.}}, {225.381, {c1 -> 1., c2 -> 0.393763, c3 -> 0., 
   c4 -> 1.}}, {312.48, {c1 -> 1., c2 -> 0., c3 -> 0., 
   c4 -> 1.}}, {298.507, {c1 -> 1., c2 -> 0., c3 -> 0., 
   c4 -> 0.483703}}, {289.22, {c1 -> 1., c2 -> 0., c3 -> 0., 
   c4 -> 0.0869823}}, {242.738, {c1 -> 0.536497, c2 -> 0., c3 -> 0., 
   c4 -> 0.}}, {206.276, {c1 -> 0., c2 -> 4.11611*10^-9, 
   c3 -> 8.75751*10^-6, c4 -> 0.}}, {143.553, {c1 -> 1., c2 -> 1., 
   c3 -> 0., c4 -> 1.}}, {196.936, {c1 -> 1., c2 -> 1., c3 -> 0., 
   c4 -> 1.}}, {245.969, {c1 -> 1., c2 -> 0.331846, c3 -> 0., 
   c4 -> 1.}}, {313.172, {c1 -> 1., c2 -> 0., c3 -> 0., 
   c4 -> 1.}}, {304.824, {c1 -> 1., c2 -> 0., c3 -> 0., 
   c4 -> 0.715149}}, {298.267, {c1 -> 1., c2 -> 0., c3 -> 0., 
   c4 -> 0.4678}}, {293.266, {c1 -> 1., c2 -> 0.0000225295, 
   c3 -> 2.52908*10^-7, c4 -> 0.262797}}, {281.43, {c1 -> 0.99607, 
   c2 -> 0.00393228, c3 -> 0.739327, 
   c4 -> 0.108152}}, {270.711, {c1 -> 0.939643, c2 -> 0.000403377, 
   c3 -> 0.999011, c4 -> 0.000087475}}, {211.763, {c1 -> 1., 
   c2 -> 0.757642, c3 -> 0., c4 -> 1.}}, {261.286, {c1 -> 1., 
   c2 -> 0.265911, c3 -> 0., c4 -> 1.}}, {313.282, {c1 -> 1., 
   c2 -> 0., c3 -> 0., c4 -> 1.}}}
| improve this answer | |
$\endgroup$
  • $\begingroup$ What is the runtime for this code segment? Also, how do I tabulate this result? (I have to include of a table of $c_{i}$s for different $t$ in my report) $\endgroup$ – Aritra Saha Jan 27 '17 at 15:05
  • $\begingroup$ This took about 10sec on my computer. For making a table I think this is relevant $\endgroup$ – gpap Jan 27 '17 at 15:08
  • $\begingroup$ I will go through the link. How do you think the increase in the number of parameters will affect the runtime? $\endgroup$ – Aritra Saha Jan 27 '17 at 15:12
  • 1
    $\begingroup$ It will increase it :) By how much depends on how many more parameters, how long a time etc. $\endgroup$ – gpap Jan 27 '17 at 15:15
  • 1
    $\begingroup$ I think that is a separate question in itself but you are essentially solving a nonlinear ode system $10^4$ times to obtain the solutions for each parameter $\endgroup$ – gpap Jan 27 '17 at 15:32

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.