4
$\begingroup$

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.

$\endgroup$
  • $\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 '15 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 '15 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 '15 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$ – cyrille.piatecki Jun 4 '16 at 22:18
  • $\begingroup$ Sounds like this could be formulated as a DAE, with the gradient used in the associated algebraic equations. $\endgroup$ – J. M. will be back soon Jan 31 '17 at 5:45
1
$\begingroup$

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


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

Mathematica graphics

$\endgroup$
0
$\begingroup$

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?

Update

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

$\endgroup$
  • $\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 '15 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.