MaplePrimes Questions

Hey there guys, was wondering if I could get some help with this - I'm pretty new to maple. 

I'm trying to take a collection of coupled ODEs I have and get maple put them in state-space form for me. I've tried a few different approaches and messed around with a few different commands but havent managed to quite make it work. At the moment I'm trying to make a DiffEquation object and then use StateSpace to get the state space representation, but I keep getting this "diff-eq is not a polynomial" line as a warning when I run the DiffEquation command and as an error when I try to conver to state space.

Could anyone tell me why I'm getting this and what the best/correct way to go about this is?

Cheers

Triple Cart derivation

with(DynamicSystems)

 

cart1eom := `m__1 `*(diff(v__1(t), t, t))+c__1*(diff(v__1(t), t))+k*(v__1(t)-v__2(t)) = F

`m__1 `*(diff(diff(v__1(t), t), t))+c__1*(diff(v__1(t), t))+k*(v__1(t)-v__2(t)) = F

(1)

cart2eom := m__2*(diff(v__2(t), t, t))+c__2*(diff(v__2(t), t))+k(v__2(t)-v__3(t)) = k*(v__1-v__2(t))

m__2*(diff(diff(v__2(t), t), t))+c__2*(diff(v__2(t), t))+k(v__2(t)-v__3(t)) = k*(v__1-v__2(t))

(2)

cart3eom := m__3*(diff(v__3(t), t, t))+c__3*(diff(v__3(t), t)) = k*(v__2-v__3(t))

m__3*(diff(diff(v__3(t), t), t))+c__3*(diff(v__3(t), t)) = k*(v__2-v__3(t))

(3)

F := alpha*k__m*k__g*V(t)/(R*r)-k__m^2*k__g^2*(diff(v__1(t), t))/(R*r^2)``

expand(cart1eom)

`m__1 `*(diff(diff(v__1(t), t), t))+c__1*(diff(v__1(t), t))+k*v__1(t)-k*v__2(t) = alpha*k__m*k__g*V(t)/(R*r)-k__m^2*k__g^2*(diff(v__1(t), t))/(R*r^2)

(4)

"(->)"

`m__1 `*(diff(diff(v__1(t), t), t))+c__1*(diff(v__1(t), t))+k*v__1(t)-k*v__2(t)-alpha*k__m*k__g*V(t)/(R*r)+k__m^2*k__g^2*(diff(v__1(t), t))/(R*r^2) = 0

(5)

collect(`m__1 `*(diff(diff(v__1(t), t), t))+c__1*(diff(v__1(t), t))+k*v__1(t)-k*v__2(t)-alpha*k__m*k__g*V(t)/(R*r)+k__m^2*k__g^2*(diff(v__1(t), t))/(R*r^2) = 0, [v__1(t), v__2, diff(v__1(t), t, t), diff(v__1(t), t)])

k*v__1(t)-k*v__2(t)+`m__1 `*(diff(diff(v__1(t), t), t))+(c__1+k__m^2*k__g^2/(R*r^2))*(diff(v__1(t), t))-alpha*k__m*k__g*V(t)/(R*r) = 0

(6)

"(->)"

cart1eom``

(7)

eoms := [cart1eom, cart2eom, cart3eom]

[`m__1 `*(diff(diff(v__1(t), t), t))+c__1*(diff(v__1(t), t))+k*(v__1(t)-v__2(t)) = alpha*k__m*k__g*V(t)/(R*r)-k__m^2*k__g^2*(diff(v__1(t), t))/(R*r^2), m__2*(diff(diff(v__2(t), t), t))+c__2*(diff(v__2(t), t))+k(v__2(t)-v__3(t)) = k*(v__1-v__2(t)), m__3*(diff(diff(v__3(t), t), t))+c__3*(diff(v__3(t), t)) = k*(v__2-v__3(t))]

(8)

NULL

sys_1 := DiffEquation(eoms, [V], [v__3], statevariable = [v__1, v__2, v__3])

_m2358947700544

(9)

StateSpace(sys_1)

_m2358947252000

(10)

NULL

Download 403TripleCart.mw

I`m trying execute the example of Deeplearnig:

with(DeepLearning);
v1 := Vector(8, i -> i, datatype = float[8]);
v2 := Vector(8, [-1.0, 1.0, 5.0, 11.0, 19.0, 29.0, 41.0, 55.0], datatype = float[8]);
model := Sequential([DenseLayer(1, inputshape = [1])]);
model := Vector(2, {(1) = Typesetting:-mi("`DeepLearning 

   Model`"), (2) = Typesetting:-mi("`<keras.engine.sequential.Se\

  quential object at 0x000001C5B6520700>`")})


model:-Compile(optimizer = "sgd", loss = "mean_squared_error");
model:-Fit(v1, v2, epochs = 500);
"<Python object: <keras.callbacks.History object at 0x000001C5CC\

  EE9DE0>>"


convert("<Python object: <keras.callbacks.History object at 0x000001C5CCEE9DE0>>", 'symbol');
<Python object: <keras.callbacks.History object at 0x000001C5CCE\

  E9DE0>>


model:-Predict([10]);

 

But, finally, there is this error:

Error, (in Predict) AttributeError: 'CatchOutErr' object has no attribute 'flush'

['Traceback (most recent call last):\n', '  File "C:\\Program Files\\Maple 2023\\Python.X86_64_WINDOWS\\lib\\site-packages\\keras\\utils\\traceback_utils.py", line 70, in error_handler\n    raise e.with_traceback(filtered_tb) from None\n', '  File "C:\\Program Files\\Maple 2023\\Python.X86_64_WINDOWS\\lib\\site-packages\\keras\\utils\\io_utils.py", line 80, in print_msg\n    sys.stdout.flush()\n', "AttributeError: 'CatchOutErr' object has no attribute 'flush'\n"]

What`s is happening?

Thanks!

I have a complicated bivariate function f(Gamma,rho) that is a RootOf of a quartic. I know that it is strictly positive (one of the four roots at least) for Gamma=0..10 and rho in (-1,+1), with bounds excluded.

I need to find the signs of its first and second derivatives (wrt to Gamma and wrt to rho: 4 derivatives in total).

I encounter numerical issues when I plot3d the derivatives using D[]() vs. fdiff() (numerical function evaluations of the RootOf). I was hoping for the two commands to produce the same output, but they don't it seems. What's going on?

Script:

restart;
_quartic := RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2);
convert(_quartic,radical):
f(Gamma,rho) := simplify(%):

RootOf((8*rho^3+24*rho^2+24*rho+8)*_Z^4+(-12*Gamma*rho^3-12*Gamma*rho^2+12*Gamma*rho+12*Gamma)*_Z^3+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*_Z^2+(4*Gamma*rho^2-4*Gamma)*_Z-Gamma^2*rho^2+2*rho*Gamma^2-Gamma^2)

(1)


Synthetic representation of derivatives

der1_Gamma := diff(_quartic, Gamma):
der1_rho := diff(_quartic, rho):

Diff('f(Gamma,rho)', Gamma) = collect~(normal(eval(der1_Gamma, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');
Diff('f(Gamma,rho)', rho) = collect~(normal(eval(der1_rho, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');

der2_Gamma := diff(der1_Gamma, Gamma):
der2_rho := diff(der1_rho, rho):

Diff('f(Gamma,rho)', Gamma$2) = collect~(normal(eval(der2_Gamma, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');
Diff('f(Gamma,rho)', rho$2) = collect~(normal(eval(der2_rho, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');

Diff(f(Gamma, rho), Gamma) = -((-6*rho^3-6*rho^2+6*rho+6)*f(Gamma, rho)^3+(5*Gamma*rho^3-5*Gamma*rho^2-5*Gamma*rho+5*Gamma)*f(Gamma, rho)^2+(2*rho^2-2)*f(Gamma, rho)-Gamma*rho^2+2*Gamma*rho-Gamma)/((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)

 

Diff(f(Gamma, rho), rho) = -(1/2)*((24*rho^2+48*rho+24)*f(Gamma, rho)^4+(-36*Gamma*rho^2-24*Gamma*rho+12*Gamma)*f(Gamma, rho)^3+(15*Gamma^2*rho^2-10*Gamma^2*rho-5*Gamma^2-8*rho-8)*f(Gamma, rho)^2+8*Gamma*rho*f(Gamma, rho)-2*rho*Gamma^2+2*Gamma^2)/((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)

 

Diff(f(Gamma, rho), Gamma, Gamma) = ((448*rho^8+1792*rho^7+1792*rho^6-1792*rho^5-4480*rho^4-1792*rho^3+1792*rho^2+1792*rho+448)*f(Gamma, rho)^8+(-1632*Gamma*rho^8-3264*Gamma*rho^7+3264*Gamma*rho^6+9792*Gamma*rho^5-9792*Gamma*rho^3-3264*Gamma*rho^2+3264*Gamma*rho+1632*Gamma)*f(Gamma, rho)^7+(2120*Gamma^2*rho^8-8480*Gamma^2*rho^6-208*rho^7+12720*Gamma^2*rho^4-624*rho^6-208*rho^5-8480*Gamma^2*rho^2+1040*rho^4+1040*rho^3+2120*Gamma^2-208*rho^2-624*rho-208)*f(Gamma, rho)^6+(-1200*Gamma^3*rho^8+2400*Gamma^3*rho^7+2400*Gamma^3*rho^6-7200*Gamma^3*rho^5+640*Gamma*rho^7+640*Gamma*rho^6+7200*Gamma^3*rho^3-1920*Gamma*rho^5-2400*Gamma^3*rho^2-1920*Gamma*rho^4-2400*Gamma^3*rho+1920*Gamma*rho^3+1200*Gamma^3+1920*Gamma*rho^2-640*Gamma*rho-640*Gamma)*f(Gamma, rho)^5+(250*Gamma^4*rho^8-1000*Gamma^4*rho^7+1000*Gamma^4*rho^6+1000*Gamma^4*rho^5-632*Gamma^2*rho^7-2500*Gamma^4*rho^4+632*Gamma^2*rho^6+1000*Gamma^4*rho^3+1896*Gamma^2*rho^5+1000*Gamma^4*rho^2-1896*Gamma^2*rho^4+16*rho^6-1000*Gamma^4*rho-1896*Gamma^2*rho^3+32*rho^5+250*Gamma^4+1896*Gamma^2*rho^2-16*rho^4+632*Gamma^2*rho-64*rho^3-632*Gamma^2-16*rho^2+32*rho+16)*f(Gamma, rho)^4+(240*Gamma^3*rho^7-720*Gamma^3*rho^6+240*Gamma^3*rho^5+1200*Gamma^3*rho^4-32*Gamma*rho^6-1200*Gamma^3*rho^3-240*Gamma^3*rho^2+96*Gamma*rho^4+720*Gamma^3*rho-240*Gamma^3-96*Gamma*rho^2+32*Gamma)*f(Gamma, rho)^3+(-25*Gamma^4*rho^7+125*Gamma^4*rho^6-225*Gamma^4*rho^5+125*Gamma^4*rho^4+125*Gamma^4*rho^3-225*Gamma^4*rho^2+125*Gamma^4*rho-25*Gamma^4)*f(Gamma, rho)^2+(16*Gamma^3*rho^6-64*Gamma^3*rho^5+80*Gamma^3*rho^4-80*Gamma^3*rho^2+64*Gamma^3*rho-16*Gamma^3)*f(Gamma, rho)-5*Gamma^4*rho^6+30*Gamma^4*rho^5-75*Gamma^4*rho^4+100*Gamma^4*rho^3-75*Gamma^4*rho^2+30*Gamma^4*rho-5*Gamma^4)/(((16*rho^2+32*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^2+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^2-10*Gamma^2*rho+5*Gamma^2-4*rho-4)*f(Gamma, rho)+2*Gamma*rho-2*Gamma)*((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)^2)

 

Diff(f(Gamma, rho), rho, rho) = (1/4)*((21504*rho^6+129024*rho^5+322560*rho^4+430080*rho^3+322560*rho^2+129024*rho+21504)*f(Gamma, rho)^10+(-80640*Gamma*rho^6-347136*Gamma*rho^5-526080*Gamma*rho^4-245760*Gamma*rho^3+157440*Gamma*rho^2+199680*Gamma*rho+56064*Gamma)*f(Gamma, rho)^9+(127680*Gamma^2*rho^6+336512*Gamma^2*rho^5+122944*Gamma^2*rho^4-319744*Gamma^2*rho^3-18176*rho^5-246976*Gamma^2*rho^2-90880*rho^4+40576*Gamma^2*rho-181760*rho^3+53696*Gamma^2-181760*rho^2-90880*rho-18176)*f(Gamma, rho)^8+(-110016*Gamma^3*rho^6-112128*Gamma^3*rho^5+191808*Gamma^3*rho^4+172032*Gamma^3*rho^3+57344*Gamma*rho^5-105792*Gamma^3*rho^2+191488*Gamma*rho^4-59904*Gamma^3*rho+192512*Gamma*rho^3+24000*Gamma^3+2048*Gamma*rho^2-94208*Gamma*rho-37888*Gamma)*f(Gamma, rho)^7+(54960*Gamma^4*rho^6-28480*Gamma^4*rho^5-102480*Gamma^4*rho^4+56960*Gamma^4*rho^3-74176*Gamma^2*rho^5+40080*Gamma^4*rho^2-126144*Gamma^2*rho^4-28480*Gamma^4*rho+41088*Gamma^2*rho^3+7440*Gamma^4+138368*Gamma^2*rho^2+5120*rho^4+19776*Gamma^2*rho+20480*rho^3-25536*Gamma^2+30720*rho^2+20480*rho+5120)*f(Gamma, rho)^6+(-15300*Gamma^5*rho^6+30000*Gamma^5*rho^5-2100*Gamma^5*rho^4-24000*Gamma^5*rho^3+50048*Gamma^3*rho^5+14100*Gamma^5*rho^2+5632*Gamma^3*rho^4-6000*Gamma^5*rho-89856*Gamma^3*rho^3+3300*Gamma^5-1024*Gamma^3*rho^2-13056*Gamma*rho^4+39808*Gamma^3*rho-31232*Gamma*rho^3-4608*Gamma^3-15360*Gamma*rho^2+10752*Gamma*rho+7936*Gamma)*f(Gamma, rho)^5+(1875*Gamma^6*rho^6-6250*Gamma^6*rho^5+7125*Gamma^6*rho^4-3500*Gamma^6*rho^3-18316*Gamma^4*rho^5+2125*Gamma^6*rho^2+25540*Gamma^4*rho^4-2250*Gamma^6*rho+12072*Gamma^4*rho^3+875*Gamma^6-26520*Gamma^4*rho^2+12992*Gamma^2*rho^4+6244*Gamma^4*rho+10240*Gamma^2*rho^3+980*Gamma^4-15488*Gamma^2*rho^2-9728*Gamma^2*rho-512*rho^3+3008*Gamma^2-1536*rho^2-1536*rho-512)*f(Gamma, rho)^4+(3320*Gamma^5*rho^5-9240*Gamma^5*rho^4+7600*Gamma^5*rho^3-560*Gamma^5*rho^2-6144*Gamma^3*rho^4-1320*Gamma^5*rho+4992*Gamma^3*rho^3+200*Gamma^5+6784*Gamma^3*rho^2-4992*Gamma^3*rho+1024*Gamma*rho^3-640*Gamma^3+1536*Gamma*rho^2-512*Gamma)*f(Gamma, rho)^3+(-200*Gamma^6*rho^5+800*Gamma^6*rho^4-1200*Gamma^6*rho^3+800*Gamma^6*rho^2+1248*Gamma^4*rho^4-200*Gamma^6*rho-3136*Gamma^4*rho^3+1920*Gamma^4*rho^2+576*Gamma^4*rho-768*Gamma^2*rho^3-608*Gamma^4+768*Gamma^2*rho)*f(Gamma, rho)^2+(-16*Gamma^5*rho^4+192*Gamma^5*rho^3-480*Gamma^5*rho^2+448*Gamma^5*rho+256*Gamma^3*rho^3-144*Gamma^5-384*Gamma^3*rho^2+128*Gamma^3)*f(Gamma, rho)-20*Gamma^6*rho^4+80*Gamma^6*rho^3-120*Gamma^6*rho^2+80*Gamma^6*rho-32*Gamma^4*rho^3-20*Gamma^6+96*Gamma^4*rho^2-96*Gamma^4*rho+32*Gamma^4)/(((16*rho^2+32*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^2+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^2-10*Gamma^2*rho+5*Gamma^2-4*rho-4)*f(Gamma, rho)+2*Gamma*rho-2*Gamma)*((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)^2)

(2)


Signs of derivatives: fdiff (numerical function evaluations of the RootOf) vs. D[]()

restart;
with(plots):

_quartic := RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2):

plot3d(_quartic, Gamma=0..10, rho=-1..+1, labels=[Gamma,rho,Lambda(Gamma,rho)],axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30]);
 

 

Define it as a f and test it for Gamma=1 and rho=0.5

f := (Gamma,rho) -> RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2):
evalf(f(1.0,0.5));

HFloat(0.5110796212870378)

(3)

Value at zero:

f(0,0):
allvalues(%):
fl := select(is, [allvalues(f(0,0))], positive)[];evalf(%);

(1/2)*2^(1/2)

 

.7071067810

(4)

Value at infinity (commented out because too slow)

#limit(f(x,y), {x = infinity, y = 0}):
#fh := select(is, [allvalues(%)], positive)[];evalf(%);

Derivative at zero:

allvalues([D[1](f)(0,0)]):
Dfl := %[1][];

-1/4

(5)

Derivative at a point, evaluated, vs numerical derivative at a point:

D[1](f)(1,0.5):
evalf(%);
fdiff(f(x,y), x, {x = 1.0, y = 0.5});
fdiff(f, [1], [1.0,0.5]);

D[2](f)(1,0.5):
evalf(%);
fdiff(f(x,y), y, {x = 1.0, y = 0.5});
fdiff(f, [2], [1.0,0.5]);

HFloat(-0.05086932918910799)

 

-0.5086932919e-1

 

-0.5086932919e-1

 

HFloat(-0.05166477232109392)

 

-0.5166477232e-1

 

-0.5166477232e-1

(6)

Can make a function out of fdiff

fDfG := (Gamma,rho) -> fdiff(f, [1], [Gamma,rho]);
fDfr := (Gamma,rho) -> fdiff(f, [2], [Gamma,rho]);

proc (Gamma, rho) options operator, arrow; fdiff(f, [1], [Gamma, rho]) end proc

 

proc (Gamma, rho) options operator, arrow; fdiff(f, [2], [Gamma, rho]) end proc

(7)

Check for numerical values close to thresholds:

Digits := 15:
evalf('D[1]'(f)(0.1e-8,0.5));fdiff(f, [1], [0.1e-8,0.5]);
evalf('D[1]'(f)(0.1e-7,0.5));fdiff(f, [1], [0.1e-7,0.5]);
evalf('D[1]'(f)(0.1e-5,0.5));fdiff(f, [1], [0.1e-5,0.5]);
evalf('D[1]'(f)(0.00001,0.5));fdiff(f, [1], [0.00001,0.5]);
evalf('D[1]'(f)(0.001,0.5));fdiff(f, [1], [0.001,0.5]);


evalf('D[2]'(f)(1,-0.99));fdiff(f, [2], [1,-0.99]);
evalf('D[2]'(f)(1,-0.97));fdiff(f, [2], [1,-0.97]);
evalf('D[2]'(f)(1,-0.1));fdiff(f, [2], [1,-0.1]);
evalf('D[2]'(f)(1,0.98));fdiff(f, [2], [1,0.98]);
evalf('D[2]'(f)(1,-0.99));fdiff(f, [2], [1,-0.99]);

57735026.8022959

 

57735026.8022959

 

-0.833333329724894e-1

 

-0.833333329724894e-1

 

-0.833332972489415e-1

 

-0.833332972489415e-1

 

-0.833329724894151e-1

 

-0.833329724894151e-1

 

-0.832972489466445e-1

 

-0.832972489466445e-1

 

-223.615892086941

 

-223.615892086941

 

-43.0236130145893

 

-43.0236130145893

 

-.212392503268663

 

-.212392503268663

 

-0.127828146340716e-2

 

-0.127828146340716e-2

 

-223.615892086941

 

-223.615892086941

(8)

Compare with D (vertical range here to prevent effect of large values from fdiff near zero):

d1G := plot3d([D[1](f), fDfG], 0..10, -0.95..+0.95, view=-0.3..0, color = [red, blue]);
d1r := plot3d([D[2](f), fDfr], 0..10, -0.95..+0.95, color = [red, blue]);

 

 

 

Second derivatives:

evalf('D[1,1]'(f)(1.0,0.5));
fdiff(f, [1, 1], [1.0,0.5]);

evalf('D[2,2]'(f)(1.0,0.5));
fdiff(f, [2, 2], [1.0,0.5]);

fD2fG := (Gamma,rho) -> fdiff(f, [1, 1], [Gamma]);
fD2fr := (Gamma,rho) -> fdiff(f, [2, 2], [Gamma]);

0.266607527050519e-1

 

0.266607527050519e-1

 

.151600577769391

 

.151600577769391

 

proc (Gamma, rho) options operator, arrow; fdiff(f, [1, 1], [Gamma]) end proc

 

proc (Gamma, rho) options operator, arrow; fdiff(f, [2, 2], [Gamma]) end proc

(9)

d2G:= plot3d([D[1,1](f), fD2fG], 0..10, -0.9..+0.9, color = [red, blue]);
d2r:= plot3d([D[2,2](f), fD2fr], 0..10, -0.9..+0.9, color = [red, blue]);
 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

d1d2G := plot3d([fDfG, fD2fG], 0.1e-6 .. 10, -0.98 .. +0.98, axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30], size=[1000,1000]);
d1d2r := plot3d([fDfr, fD2fr], 0.1e-6 .. 10, -0.98 .. +0.98, axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30], size=[1000,1000]);

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 
 

NULL

Download signs_derivatves_bivariate.mw

Why does no substitution work on functions with s(n+1,t) (see eqns (2-8))? Also (&PartialD;)/(&PartialD; (sigma2*t))=1/(sigma2)(&PartialD;)/(&PartialD; t), how can I do it on Maple?shift.mw

I am trying to load the third-party package 'CPC Program Library (qub.ac.uk)' by following the instructions as in 'how to install wkptest? - MaplePrimes'. But encounter with Errors: '

Error, `:` unexpected
with(wkptest);
Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received wkptest

restart:

sourcefolder:=cat(kernelopts('C:\Users\ahmed\Downloads\adty_v1_0'),"/wkptest");
installfolder:=cat(kernelopts('homedir'),"/maple/toolbox/wkptest/lib");
FileTools:-MakeDirectory(installfolder, 'recurse'=true);
libraryfile:=cat(installfolder,"/wkptest.mla");
try
  FileTools:-Remove(libraryfile);
catch:
end try:
LibraryTools:-Create(libraryfile);
libname:=libraryfile,libname;
read cat(sourcefolder,"/wkptest_cpc");

Error, `:` unexpected

 

with(wkptest);

Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received wkptest

 
 

 

Download exam_wkptest1.mws

When I calculate the difference between two dates with the same day using the DateDifference Calendar function, I don't get an integer number of months, but results with extra days, hours, minutes, seconds and milliseconds. How can I calulate the actual number of months?

I appreciate any help.

restart

with(Calendar)

d1 := Date(2024, 1, 5, 0, 0, 0)

_m4674925856

(1)

d2 := Date(2024, 5, 5, 0, 0, 0)

_m4881744064

(2)

DateDifference(d1, d2, 'units' = 'mixed')

4*Units:-Unit(mo)+2*Units:-Unit(d)+20*Units:-Unit(h)+3*Units:-Unit(min)+48*Units:-Unit(s)+800*Units:-Unit(ms)

(3)

NULL

NULL

 

NULL

Download DateDifference.mw

How can I calculate the integral of Legendre function P(n,theta) ?

int(cos(theta)*P(n,theta),theta=0..pi) ;

I never used Maple's StringTools:-RegSubs before.

I have a Latex string generated by Maple that has this form

"......    \\textrm{ ......   }  .... "

Where the dots mean anything, including space characters.  I need to change the part of the string  "\\textrm{ ...... }" to this

"......  \\begin{minipage}{\\linewidth}\\textrm{  ......   }\\end{minipage}  .... "

It is easy to change the opening part. Here is an example

s:="A& B& \\textrm{The fundamental matrix has } &C ";
s:=StringTools:-SubstituteAll(s,"\\textrm{","\\begin{minipage}{\\linewidth}\\textrm{");

Gives

The problem is how to replace the closing "}"  with "}\\end{minipage}".   I can't just replace "}" with "}\\end{minipage}" using SubstituteAll because I need to only change the "}" that closes the corresponding "\\textrm{" and not any other "}" that could be in the string.

So I need to use pattern matching or regular expression substitutions. The examples in help are not easy to understand.

 

Basically I need regular expression that matches 

                        "\\textrm{ WILDCARD * }" 

and change it to 

                        "\\begin{minipage}{\\linewidth}\\textrm{ WILDCARD * }\\end{minipage}" 

Any one know how to use Maple's egular expression to do this substitution using the above example?

Maple 2024 on windows 10.

Update

Thanks for the replies. I ended up writing small function that simply scan the string and do the replacement. 

Hi,

I want to analyze a quartic equation: rho-analysis.mw. I am interested in positive roots, which I need expressed in explicit/closed form.

I include four questions in the script.

Thanks!

EDIT (since my question was originally tagged as incomplete and even duplicate):
Above all I am trying to understand:

  1. Why implicitplot3d returns a unique strictly positive root for all rho and Gamma but plot() of the quartic evaluated for a specific pair of value returns three positive roots?
  2. rho is a correlation coefficient, thus bounded in (-1,+1): how to adjust the precision of the slider in Explore() so that I can play around with multiple rho values within said range?

Can this be better done in Maple ? , see worksheet.

Printlevel can be used to understand how Maple processes input. In some cases printlevel has to be increased to high values which produces an ennormeous amount of output.

In case a particular procedure has been identified for detailed analysis it would be desirable only for this procedure to set a printlevel. Kind of a local printlevel.

Does something like this or other ways exist?

trace as I understand it only traces the procedure given as argument but not levels "below".

In this exercise, we are going back to the simple SIR model, without births or deaths, to look at the effect of vaccination. The aim of this activity is to represent vaccination in a very simple way - we are assuming it already happened before we run our model! By changing the initial conditions, we can prepare the population so that it has received a certain coverage of vaccination.

We are starting with the transmission and recovery parameters b = 0.4 days-1 and  c=0.1 days-1. To incorporate immunity from vaccination in the model, we assume that a proportion p of the total population starts in the recovered compartment, representing the vaccine coverage and assuming the vaccine is perfectly effective. Again, we assume the epidemic starts with a single infected case introduced into the population.​
We are going to model this scenario for a duration of 2 years, assuming that the vaccine coverage is 50%, and plot the prevalence in each compartment over time.

I have used the following code to solve this problem but I face difficulties when I am trying to define and plot the effective reproduction number.

Any suggestions? Thank you in advance!
with(Physics, diff);
with(LinearAlgebra);

with(plots);
with(DEtools);

b := 0.4;
c := 0.1;
n := 10^6;
p := 0.5;
deS := diff(S(t), t) = -b*S(t)*I0(t);
deI := diff(I0(t), t) = b*S(t)*I0(t) - c*I0(t);
deR := diff(R(t), t) = c*I0(t);

F := dsolve([deS, deI, deR, S(0) = 1 - p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "Proportion of Population"], title = "SIR Model with vaccination")


#define the effective reproduction number
Reff := t -> b*1/c*S(t)*1/n;
Reff(100);

I have been trying to solve the following diffusion equations

  diff(u(z, t), t) = D*diff(u(z, t), z, z) + A*z*exp(lambda*z)

my boundary canditions are

  u(0,t)=0,D(u)(H,t)=0,u(initinity,t)=0

I have followed the examples on tte link below, but I am unable to get the solution when I write the boundary conditions. It there anything I am doing wrongly? I would be most grateful if an email can be sent to me with with the solution that I can download fram an attachment

Thanks

https://www*mapleprimes*com/posts/209970-Exact-Solutions-For-PDE-And-Boundary--Initial-Conditions-2018

SI_MODEL.mw
I have tried the following code to simulate the following SI Model. But when I run it, it gives me several mistakes and most importantly does not calculate the Jacobian Matrix and the eigenvalues. 

Any help would be deeply appreciated


Many of the most interesting dynamics in nature have to do with interactions between organisms. These interactions are often subtle, indirect, and difficult to detect. Interactions in which one organism consumes all or part of another. This includes predator-prey, herbivore-plant, and parasite-host interactions. These linkages are the prime movers of energy through food chains. They are an important factor in the ecology of populations, determining the mortality of prey and the birth of new predators.

Mathematical models and logic suggest that a coupled system of predator and prey should cycle: predators increase when prey is abundant, prey are driven to low numbers by predation, the predators decline, and the prey recover, ad infinitum. One such model that simulates predator-prey interactions is the Lotka - Volterra Model.

 


We will discuss a behavior of 2 D - subsystems: SI - model.


SI Model (without predator)

 

We will study a mathematical model that appears in eco - eco-epidemiology.


The SI model describes the interactions between S-prey with density x and I-prey with density y under the assumption z≡0


Break the population into compartments.

• 

Prey (S)

• 

Infected prey (I)


SI Equations :

(D(x))(t) = r*x(t)-b*x(t)^2-c*x(t)*y(t)-beta*x(t)*y[t]/(a+x(t))

(D(y))(t) = -mu*y(t)+beta*x(t)*y(t)/(a+x(t))

– 

 r stands for the intrinsic growth rate of S-prey.

– 

b defines the intra-class competition in S-prey.

– 

c characterizes the inter-class competition between S-prey and I-prey.

– 

μ stands for the mortality of the I-prey

– 

β considered as the bifurcation parameter.

• 

Note: The dynamics of the system depend on parameter β

 

Modeling with Differential Equations

 

• 

We will fix the parameters in the study to have the following values r = 1, b = 1, c = 1/100, μ = 4/10, a=1/2.75.

• 

The initial values will be x(0)=0.2, y(0)=0.05


As mentioned above the dynamics of the system depend on parameter β. So we will consider two cases when β=0.75 and β=1


We will solve the system for the above values of parameters and initial conditions and for a time interval [0,200] , then we will create its plots

 

Case 1

  restart

beta := .75; r := 1; b := 1; c := 1/100; mu := 4*(1/10); a := 1/2.75; f := r*x(t)-b*x(t)^2-c*x(t)*y(t)-(3/4)*x(t)*y(t)/(a+x(t)); g := -mu*y(t)+(3/4)*x(t)*y(t)/(a+x(t)); deq1 := diff(x(t), t) = f; deq2 := diff(y(t), t) = g; equilibrio := solve({f = 0, g = 0}, {x(t), y(t)}); jacobian := simplify(Matrix(2, 2, [[diff(f, x(t)), diff(f, y(t))], [diff(g, x(t)), diff(g, y(t))]])); A := simplify(eval(jacobian, equilibrio[1])); B := simplify(eval(jacobian, equilibrio[2])); Q := simplify(eval(jacobian, equilibrio[3])); eigen1 := Eigenvalues(A); eigen2 := Eigenvalues(B); eigen3 := Eigenvalues(Q); soln := dsolve({deq1, deq2, x(0) = .2, y(0) = 0.5e-1}, numeric, output = listprocedure); plots:-display(plots:-odeplot(soln, [t, x(t)], 0 .. 200, color = blue, legend = ["x(t)"]), plots:-odeplot(soln, [t, y(t)], 0 .. 200, color = red, legend = ["y(t)"]), labels = ["t", "Population"], title = "Population Dynamics")

 

NULL

 

 

Case 2

restart``

beta := 1; r := 1; b := 1; c := 1/100; Mu := 4*(1/10); a := 1/2.75; f := r*x(t)-b*x(t)^2-c*x(t)*y(t)-x(t)*y(t)/(a+x(t)); g := -Mu*y(t)+x(t)*y(t)/(a+x(t)); deq1 := diff(x(t), t) = f; deq2 := diff(y(t), t) = g; equilibrio := solve({f = 0, g = 0}, {x(t), y(t)}); jacobian := Matrix(2, 2, [[diff(f, x(t)), diff(f, y(t))], [diff(g, x(t)), diff(g, y(t))]], simplify); J_at_equilibrio := [eval(jacobian, equilibrio)]; eigen1 := Eigenvalues(A); eigen2 := Eigenvalues(B); eigen3 := Eigenvalues(Q); soln := dsolve({deq1, deq2, x(0) = .2, y(0) = 0.5e-1}, numeric); plots:-display(plots:-odeplot(soln, [t, x(t)], 0 .. 200, color = blue, legend = ["x(t)"]), plots:-odeplot(soln, [t, y(t)], 0 .. 200, color = red, legend = ["y(t)"]), labels = ["t", "Population"], title = "Population Dynamics with Beta = 1")

 

NULL


 

Download SI_MODEL.mw

 

 

First 109 110 111 112 113 114 115 Last Page 111 of 2428