I am trying to compute the maximum value of the solution of a system with two ODEs using Maple. I have firstly solved the system itself:
> with(DEtools):with(plots):
> a1:=0.00875;a2:=0.075;b1:=7.5;b2:=2.5;d1:=0.0001;d2:=0.0001;g:=4*10^(-8);K1:=5000;K2:=2500;n:=2;m:=2;
> dsol:= dsolve({
diff( x(t), t ) = a1+b1*x(t)^n/(K1^n+x(t)^n)-g*x(t)*y(t)-d1*x(t),
diff( y(t), t ) = a2+b2*x(t)^m/(K2^m+x(t)^m)-d2*y(t),
x(0) = 1000, y(0) = 1000}, numeric, output = listprocedure);
> xt:= eval( x(t), dsol );
yt:= eval( y(t), dsol );
> X:=plot(xt(t),t=0..50000,color=blue,legend="x(t)"):
Y:=plot(yt(t),t=0..50000,color=green,legend="y(t)"):
> display([X,Y]);
I got the solutions of the system on xt and yt, but they are numeric solutions. Therefore, the function of Maple maximize() does not work:
> maximize(xt);
> maximize(xt(t),t=0..20000);
Is it possible to compute the maximum value of a numeric function with Maple?
Your two curves xt
and yt
each have a single local maximum within your range of t=0..50000
, so you can use the Optimization
package for this in a straightforward way.
restart;
with(plots):
a1:=0.00875: a2:=0.075: b1:=7.5: b2:=2.5: d1:=0.0001:
d2:=0.0001: g:=4*10^(-8): K1:=5000: K2:=2500: n:=2: m:=2:
dsol:= dsolve({diff(x(t),t)=a1+b1*x(t)^n/(K1^n+x(t)^n)-g*x(t)*y(t)-d1*x(t),
diff(y(t),t)=a2+b2*x(t)^m/(K2^m+x(t)^m)-d2*y(t),
x(0)=1000, y(0)=1000}, numeric, output=listprocedure):
xt:= eval(x(t), dsol):
yt:= eval(y(t), dsol):
X:=plot(xt(t), t=0..50000, color=blue, legend="x(t)"):
Y:=plot(yt(t), t=0..50000, color=green, legend="y(t)"):
xmax:=Optimization:-Maximize(xt, 0..50000):
[xmax[2][1],xmax[1]];
[9460.78688552799, 11193.0618953179]
ymax:=Optimization:-Maximize(yt, 0..50000):
[ymax[2][1],ymax[1]];
[21471.8648785947, 19006.6009784691]
display( Y, pointplot([[ymax[2][1],ymax[1]]], symbolsize=20),
X, pointplot([[xmax[2][1],xmax[1]]], symbolsize=20) );
For your simple example this works out pretty well.
If your xt
or yt
had many local maxima you could try calling Maximize
with its method=branchandbound
option.
And then there are other approaches where you could augment your DE system with new dependent variables like, say, xd(t)=diff(x(t),t)
(along with appropriate ICs) and get dsolve/numeric itself to notice when that becomes zero (an extreme point) using the events
facilities, or use fsolve
on it.