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.

  • $\begingroup$ I guess I cannot help you unless you specify what x' mean? D[x[t],t]? D[x[a],a]? ?? $\endgroup$
    – chris
    Aug 9, 2015 at 12:49
  • 1
    $\begingroup$ 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? $\endgroup$
    – bbgodfrey
    Aug 10, 2015 at 21:38
  • $\begingroup$ 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). $\endgroup$
    – bbgodfrey
    Aug 10, 2015 at 21:51
  • $\begingroup$ 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 $\endgroup$ Jun 4, 2016 at 22:18
  • $\begingroup$ Sounds like this could be formulated as a DAE, with the gradient used in the associated algebraic equations. $\endgroup$ Jan 31, 2017 at 5:45

2 Answers 2


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

Mathematica graphics

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]

Mathematica graphics


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

Mathematica graphics


This is not an answer but a comment

How about

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

Mathematica graphics

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


oops I had not seen that your optimization problem involves 'x,y'…

  • $\begingroup$ well, maximize 3a^0.8b^0.2+a^2/(a-1)+b^2/(b-1) in 0<=x<=0.99999, 0<=y<=0.99999 returns in a compact set a=0.43917, b=0.27115. Further, we get derivatives x', y', and attain a new state of the system, and once again calculate a, b via optimization. since for every x,y there is a unique solution (a,b), global extreme, it may be said a=f(x,y),b=g(x,y). here a question arises, do we know derivatives? - yes, do we have initial values? - yes. thus there should be an integral curve. $\endgroup$
    – adam
    Aug 9, 2015 at 15:47

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.