2016-03-19 14 views
2

Maple kullanarak iki ODE ile bir sistemin çözümünün maksimum değerini hesaplamaya çalışıyorum. Öncelikle sistem kendisi çözdük:Maple, sayısal bir fonksiyonun maksimum değeri

> 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]); 

Ben xt ve yt üzerine sistemin çözüm var, ama sayısal çözümleridir. Bu nedenle, Akçaağaç üst düzeye çıkarma function() çalışmaz:

> maximize(xt); 
> maximize(xt(t),t=0..20000); 

o Maple ile sayısal fonksiyonun maksimum değerini hesaplamak mümkün mü?

cevap

2

Sizin iki eğri xt ve yt her t=0..50000 aralığınız içinde tek yerel maksimum var, yani basit bir şekilde bunun için Optimization paketini kullanabilirsiniz.

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)); 
sizin basit Örneğin

enter image description here

bu oldukça iyi çalışıyor.

senin xt veya yt birçok yerel maksimum olsaydı onun method=branchandbound seçeneğiyle Maximize arama deneyebilirsiniz.

Ve sonra DE sisteminizi xd(t)=diff(x(t),t) (uygun IC'lerle birlikte) gibi yeni bağımlı değişkenler ile artırabileceğiniz ve bu sıfıra ulaştığında fark etmek için dsolve/numeric'in kendisini (başka bir nokta) kullanabileceği başka yaklaşımlar vardır. events tesislerinde veya üzerinde fsolve kullanın.

+0

Mükemmel! Çok iyi bir açıklama ve gerçekten kullanışlı kod! Çok teşekkür ederim!! –

İlgili konular