Coupled ODEs with parameters subject to optimization problem

I have an equation, but I don’t know how to program it in Mathematica:

x[0]==1; y[0]==1; x'[t]==(a+b) x[t] y[t]; y'[t]==a b y[t]/x[t];


where a, b are numerical solution of optimization problem, or system of algebraic equations:

Maximize[{3a^0.8b^0.2+a^2/(a-x)+b^2/(b-y)}, {a, b}]


If there had been an exact solution for a, b it would be possible to plug their values into the differential equation, but in this case, first the magnitudes of a and b are calculated, and only then derivatives are determined.

• I guess I cannot help you unless you specify what x' mean? D[x[t],t]? D[x[a],a]? ?? – chris Aug 9 '15 at 12:49
• Since x and y are functions of t as well as of a and b, it seems to me that Maximize must be done with respect to {a, b, t}. Is that your intent? – bbgodfrey Aug 10 '15 at 21:38
• Note also that the quantity to be Maximized becomes infinite as t approaches Log[(a^2 (1 + b))/(a + b)]/(a + b + a b) or Log[(a b^2)/(a + b - b^2)]/(a + b + a b). – bbgodfrey Aug 10 '15 at 21:51
• I am not sure but as x[0] and y[0] are given you can calculate a[x[0]] and b[x[0]]. So if you discretize your differential system of equations you can obtain x[0+h],y[0+h] and then use those two values to calculate a[x[h]], b[x[h]] and so on. In choosing a small gap h you certainly approach the solution. If you cannot manage it I will try later – cyrille.piatecki Jun 4 '16 at 22:18
• Sounds like this could be formulated as a DAE, with the gradient used in the associated algebraic equations. – J. M. will be back soon Jan 31 '17 at 5:45

First of all, I'd like to point out the function you chose for optimization doesn't have a finite maximum:

With[{k = 10^20},
Plot3D[3 a^(8/10) b^(2/10) + a^2/(a - x) + b^2/(b - y) /. x | y -> 1 //
Evaluate, {a, -k, k}, {b, -k, k}]]


Still, generally it's possible to do what you want, here I'll use

Maximize[-(a - Sin[x])^2 - (b - Sin[y])^2, {a, b}]
(* {0, {a -> Sin[x], b -> Sin[y]}} *)


as an example. As you can see this optimization problem has an symbolic solution so we can easily check if the obtained solution is correct.

Back to our problem, what we need is just a black box function:

blackbox[{x_?NumericQ, y_}] :=
Module[{a, b}, {a, b} = NArgMax[-(a - Sin[x])^2 - (b - Sin[y])^2, {a, b}]; {(a + b) x y,
a b y/x}]

solxy = NDSolveValue[{xy'[t] == blackbox[xy@t], xy[0] == {1, 1}},
xy, {t, 0, 1}]; // AbsoluteTiming
(* {3.044525, Null} *)

{solx, soly} =
ListInterpolation[#, Flatten@solxy["Grid"]] & /@ Transpose@solxy["ValuesOnGrid"];

(* Result for comparison *)
{solx2, soly2} =
NDSolveValue[{x[0] == 1, y[0] == 1, x'[t] == (a + b) x[t] y[t],
y'[t] == a b y[t]/x[t]} /.
Last@Maximize[-(a - Sin[x@t])^2 - (b - Sin[y@t])^2, {a, b}], {x, y}, {t, 0, 1}];

Plot[{solx@t, soly@t, solx2@t, soly2@t}, {t, 0, 1},
PlotStyle -> {{Thick, Dashed}, {Thick, Dashed}, Thin, Thin}, PlotRange -> All]


Update

Just saw your comment under Chris' answer. With the constraint 0 < a < 1, 0 < b < 1 there's indeed a finite maximum. Then we just need to define the black box function as:

blackbox[{x_?NumericQ, y_}] :=
Module[{a, b}, {a, b} =
NArgMax[{3 a^(8/10) b^(2/10) + a^2/(a - x) + b^2/(b - y), 0 < a < 1, 0 < b < 1}, {a,
b}]; {(a + b) x y, a b y/x}]


Rest part of the code is unchanged. The obtained result is:

Plot[{solx[t], soly[t]}, {t, 0, 1}]


This is not an answer but a comment

 {X,Y}=DSolveValue[{x[0] == 1, y[0] == 1, x'[t] == (a + b) x[t] y[t],

Then you can specify a posteriori your values for a and b?