Maple 2018 Questions and Posts

These are Posts and Questions associated with the product, Maple 2018

Someone please help me with the computation of the right and left eigenvectors. my system of equation is attached below

with(VectorCalculus):

 

interface(imaginaryunit = I)

I

(2)

I

I

(3)

sqrt(-4)

2*I

(4)

NULL

``

Limit(N(t) = N__0*exp(-mu*t)+exp(mu*t)*K/mu, t = infinity)

 

limit(N(t), t = infinity) = limit(N__0*exp(-mu*t)+exp(mu*t)*K/mu, t = infinity)

(5)

 

NULL

#to calculate the  disease free equilibrium,

NULL

E1 := -S*µ__C+`Λ__p`

-S*µ__C+Lambda__p

(6)

NULL

``

(7)

E3 := -S__A*µ__A+`Λ__A`

-S__A*µ__A+Lambda__A

(8)

NULL

``

(9)

NULL

``

(10)

NULL

solve({E1 = 0, E3 = 0}, {S, S__A})

{S = Lambda__p/µ__C, S__A = Lambda__A/µ__A}

(11)

NULL

NULL#to calculate the Endemic Equilibrium state,

Typesetting:-mparsed()

(12)

restart

with(VectorCalculus):

 

interface(imaginaryunit = I)

I

(14)

I

I

(15)

sqrt(-4)

2*I

(16)

``

E1 := `Λ__p`-(`ϕ`*`θ__B`*I__A/N__p+µ__C)*S+`ω__B`*I__B

Lambda__p-(varphi*theta__B*I__A/N__p+µ__C)*S+omega__B*I__B

(17)

E2 := `ϕ`*`θ__B`*I__A*S/N__p-`ω__B`*I__B-(`σ__B`+µ__C)*I__B

varphi*theta__B*I__A*S/N__p-omega__B*I__B-(sigma__B+µ__C)*I__B

(18)

``

(19)

E3 := `Λ__A`-(µ__A+`ϕ`*`α__B`*I__B/N__p)*S__A+`δ__A`*I__A

Lambda__A-(µ__A+varphi*alpha__B*I__B/N__p)*S__A+delta__A*I__A

(20)

E4 := `ϕ`*`α__B`*I__B*S__A/N__p-(µ__A+`δ__A`)*I__A

varphi*alpha__B*I__B*S__A/N__p-(µ__A+delta__A)*I__A

(21)

NULL

``

(22)

NULL

``

(23)

solve({E1 = 0, E2 = 0, E3 = 0, E4 = 0}, {I__A, I__B, S, S__A})

{I__A = 0, I__B = 0, S = Lambda__p/µ__C, S__A = Lambda__A/µ__A}, {I__A = -(N__p^2*µ__A^2*µ__C^2+N__p^2*µ__A^2*µ__C*omega__B+N__p^2*µ__A^2*µ__C*sigma__B+N__p^2*µ__A*µ__C^2*delta__A+N__p^2*µ__A*µ__C*delta__A*omega__B+N__p^2*µ__A*µ__C*delta__A*sigma__B-varphi^2*Lambda__A*Lambda__p*alpha__B*theta__B)/(varphi*µ__A*theta__B*(N__p*µ__A*µ__C+N__p*µ__A*sigma__B+N__p*µ__C*delta__A+N__p*delta__A*sigma__B+varphi*Lambda__p*alpha__B)), I__B = -(N__p^2*µ__A^2*µ__C^2+N__p^2*µ__A^2*µ__C*omega__B+N__p^2*µ__A^2*µ__C*sigma__B+N__p^2*µ__A*µ__C^2*delta__A+N__p^2*µ__A*µ__C*delta__A*omega__B+N__p^2*µ__A*µ__C*delta__A*sigma__B-varphi^2*Lambda__A*Lambda__p*alpha__B*theta__B)/(alpha__B*(N__p*µ__A*µ__C^2+N__p*µ__A*µ__C*omega__B+N__p*µ__A*µ__C*sigma__B+varphi*µ__C*Lambda__A*theta__B+varphi*Lambda__A*sigma__B*theta__B)*varphi), S = (N__p*µ__A*µ__C+N__p*µ__A*sigma__B+N__p*µ__C*delta__A+N__p*delta__A*sigma__B+varphi*Lambda__p*alpha__B)*µ__A*N__p*(µ__C+omega__B+sigma__B)/(alpha__B*varphi*(N__p*µ__A*µ__C^2+N__p*µ__A*µ__C*omega__B+N__p*µ__A*µ__C*sigma__B+varphi*µ__C*Lambda__A*theta__B+varphi*Lambda__A*sigma__B*theta__B)), S__A = N__p*(N__p*µ__A^2*µ__C^2+N__p*µ__A^2*µ__C*omega__B+N__p*µ__A^2*µ__C*sigma__B+N__p*µ__A*µ__C^2*delta__A+N__p*µ__A*µ__C*delta__A*omega__B+N__p*µ__A*µ__C*delta__A*sigma__B+varphi*µ__A*µ__C*Lambda__A*theta__B+varphi*µ__A*Lambda__A*sigma__B*theta__B+varphi*µ__C*Lambda__A*delta__A*theta__B+varphi*Lambda__A*delta__A*sigma__B*theta__B)/(varphi*µ__A*theta__B*(N__p*µ__A*µ__C+N__p*µ__A*sigma__B+N__p*µ__C*delta__A+N__p*delta__A*sigma__B+varphi*Lambda__p*alpha__B))}

(24)

``

J := Jacobian([E1, E2, E3, E4], [S, I__B, S__A, I__A])

Matrix(%id = 18446746854857131062)

(25)

NULL

restart

J := Matrix(4, 4, {(1, 1) = -`ϕ`*`θ__B`*I__A/N__p-µ__C, (1, 2) = `ω__B`, (1, 3) = 0, (1, 4) = -`ϕ`*`θ__B`*S/N__p, (2, 1) = `ϕ`*`θ__B`*I__A/N__p, (2, 2) = -`ω__B`-`σ__B`-µ__C, (2, 3) = 0, (2, 4) = `ϕ`*`θ__B`*S/N__p, (3, 1) = 0, (3, 2) = -`ϕ`*`α__B`*S__A/N__p, (3, 3) = -µ__A-`ϕ`*`α__B`*I__B/N__p, (3, 4) = `δ__A`, (4, 1) = 0, (4, 2) = `ϕ`*`α__B`*S__A/N__p, (4, 3) = `ϕ`*`α__B`*I__B/N__p, (4, 4) = -µ__A-`δ__A`})

Matrix(%id = 18446746579340105118)

(26)

S := `Λ__p`/µ__C

Lambda__p/µ__C

(27)

S__A := `Λ__A`/µ__A

Lambda__A/µ__A

(28)

I__B := 0

0

(29)

I__A := 0

0

(30)

NULL

0

(31)

J := Matrix(4, 4, {(1, 1) = -`ϕ`*`θ__B`*I__A/N__p-µ__C, (1, 2) = `ω__B`, (1, 3) = 0, (1, 4) = -`ϕ`*`θ__B`*S/N__p, (2, 1) = `ϕ`*`θ__B`*I__A/N__p, (2, 2) = -`ω__B`-`σ__B`-µ__C, (2, 3) = 0, (2, 4) = `ϕ`*`θ__B`*S/N__p, (3, 1) = 0, (3, 2) = -`ϕ`*`α__B`*S__A/N__p, (3, 3) = -µ__A-`ϕ`*`α__B`*I__B/N__p, (3, 4) = `δ__A`, (4, 1) = 0, (4, 2) = `ϕ`*`α__B`*S__A/N__p, (4, 3) = `ϕ`*`α__B`*I__B/N__p, (4, 4) = -µ__A-`δ__A`})

Matrix(%id = 18446746579340107518)

(32)

J := Matrix(4, 4, {(1, 1) = -µ__C, (1, 2) = `ω__B`, (1, 3) = 0, (1, 4) = -`β__1`, (2, 1) = 0, (2, 2) = -`ω__B`-`σ__B`-µ__C, (2, 3) = 0, (2, 4) = -`β__1`, (3, 1) = 0, (3, 2) = -`β__2`, (3, 3) = -µ__A, (3, 4) = `δ__A`, (4, 1) = 0, (4, 2) = `β__2`, (4, 3) = 0, (4, 4) = -µ__A-`δ__A`})

Matrix(%id = 18446746579417403630)

(33)

"simplify( ? )"

Matrix(%id = 18446746579305905318)

(34)

"LinearAlgebra:-Eigenvalues( ? )"

Vector[column](%id = 18446746579445964182)

(35)

"LinearAlgebra:-CharacteristicPolynomial( ?, lambda )"

lambda^4+(2*µ__A+delta__A+omega__B+sigma__B+2*µ__C)*lambda^3+(beta__1*beta__2+µ__A^2+4*µ__A*µ__C+µ__A*delta__A+2*µ__A*omega__B+2*µ__A*sigma__B+µ__C^2+2*µ__C*delta__A+µ__C*omega__B+µ__C*sigma__B+delta__A*omega__B+delta__A*sigma__B)*lambda^2+(beta__1*beta__2*µ__A+beta__1*beta__2*µ__C+2*µ__A^2*µ__C+µ__A^2*omega__B+µ__A^2*sigma__B+2*µ__A*µ__C^2+2*µ__A*µ__C*delta__A+2*µ__A*µ__C*omega__B+2*µ__A*µ__C*sigma__B+µ__A*delta__A*omega__B+µ__A*delta__A*sigma__B+µ__C^2*delta__A+µ__C*delta__A*omega__B+µ__C*delta__A*sigma__B)*lambda+beta__1*beta__2*µ__A*µ__C+µ__A^2*µ__C^2+µ__A^2*µ__C*omega__B+µ__A^2*µ__C*sigma__B+µ__A*µ__C^2*delta__A+µ__A*µ__C*delta__A*omega__B+µ__A*µ__C*delta__A*sigma__B

(36)

NULL

"(->)"

Vector[column](%id = 18446746579340117046)

(37)

# to find the trace we

 

Matrix(7, 7, {(1, 1) = -beta*lambda-v__1-µ, (1, 2) = v__2, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (2, 1) = v__1, (2, 2) = beta*(w-1)*lambda-µ-v__2-alpha, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (3, 1) = 0, (3, 2) = alpha, (3, 3) = -µ, (3, 4) = 0, (3, 5) = `ρ__A`, (3, 6) = `ρ__F`, (3, 7) = -(-1+k)*`ρ__Q`, (4, 1) = beta*lambda, (4, 2) = -beta*(w-1)*lambda, (4, 3) = 0, (4, 4) = -q__E-delta-µ, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = a*delta, (5, 5) = -`ρ__A`-q__A-µ, (5, 6) = 0, (5, 7) = k*`ρ__Q`, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = -delta*(-1+a), (6, 5) = 0, (6, 6) = -`ρ__F`-q__F-`δ__F`-µ, (6, 7) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = q__E, (7, 5) = q__A, (7, 6) = q__F, (7, 7) = -`ρ__Q`-`δ__Q`-µ})

Matrix(%id = 36893490965935089652)

(38)

"(->)"

-beta*lambda-v__1-7*µ+beta*(w-1)*lambda-v__2-alpha-q__E-delta-rho__A-q__A-rho__F-q__F-delta__F-rho__Q-delta__Q

(39)

 

#this shows that trace is negative

 

#to Achieve stability, the value below must be less than zero

 

(-q__E-delta-µ)*(-`&rho;__F`-q__F-`&delta;__F`-µ)*(-k*q__A*`&rho;__Q`+q__A*µ+q__A*`&delta;__Q`+q__A*`&rho;__Q`+µ^2+µ*`&delta;__Q`+µ*`&rho;__A`+µ*`&rho;__Q`+`&delta;__Q`*`&rho;__A`+`&rho;__A`*`&rho;__Q`)*µ < 0

(-q__E-delta-µ)*(-rho__F-q__F-delta__F-µ)*(-k*q__A*rho__Q+q__A*rho__Q+q__A*µ+q__A*delta__Q+rho__A*rho__Q+rho__A*µ+rho__A*delta__Q+rho__Q*µ+µ^2+µ*delta__Q)*µ < 0

(40)

 NULL

M := diff(N(t), t) = Pi-µ*N(t)

diff(N(t), t) = Pi-µ*N(t)

(41)

dsolve({M}, N(t))

{N(t) = Pi/µ+exp(-µ*t)*_C1}

(42)

eval({N(t) = Pi/µ+exp(-µ*t)*_C1}, [t = infinity])

{N(infinity) = Pi/µ+exp(-µ*infinity)*_C1}

(43)

value(%)

{N(infinity) = Pi/µ+exp(-µ*infinity)*_C1}

(44)

Limit(N(t) = Pi/µ+exp(-µ*t)*_C1, t = infinity); value(%)

Limit(N(t) = Pi/µ+exp(-µ*t)*_C1, t = infinity)

 

limit(N(t), t = infinity) = limit(Pi/µ+exp(-µ*t)*_C1, t = infinity)

(45)

 

Subs := diff(S(t), t) = -(beta*lambda+v__1+µ)*S(t)

diff(S(t), t) = -(beta*lambda+v__1+µ)*S(t)

(46)

dsolve({Subs}, S(t))

{S(t) = _C1*exp(-(beta*lambda+v__1+µ)*t)}

(47)
 

``

Download Cotton_DFE_and_Jacobian.mw

Please, I am encountering error trying to run these codes for sensitivity analysis using the formula for sensitivity analysis

``

restart;

#
# Set up numerical values for all problem parameters
#
  params:=[ Lambda__p=100,         gamma__B=0.05,      gamma__B=0.05,
                 gamma__C=0.01, omega__C=0.001,  omega__B=0.001,
            sigma__B=0.0001,     sigma__C=0.01, sigma__BC=0.01,
                theta__B=0.8,     theta__C=0.5,      mu__C=1.0,
              Lambda__A=1.0,       Lambda__w=1.0,   varphi__8.33,
            mu__A=1.0, mu__w=1.0, alpha__B=0.005, alpha__C=0.005, alpha__BC=0.15, Zeta__B=0.5, Zeta__C=0.5, delta__A=0.66, delta__w=1.33
          ]:

#
# Define main function
#
  R:= (varphi^2*theta__B*Lambda__p*alpha__B*Lambda__A)/((mu__c*mu__A*N__p^2)*(mu__A*mu__c+mu__A*omega__B+mu__A*sigma__B+mu__c*delta__A+delta__A*omega__B+delta__A*sigma__B));

varphi^2*theta__B*Lambda__p*alpha__B*Lambda__A/(mu__c*mu__A*N__p^2*(mu__A*mu__c+mu__A*omega__B+mu__A*sigma__B+mu__c*delta__A+delta__A*omega__B+delta__A*sigma__B))

(1)

#
# Compute "all" derivatives and evaluate numerically.
#
# For the purposes of this calculation "all"
# derivatives, means the derivatives with respect to
# every variable returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( diff( R, varName), params )
# ]
#
 [ seq( [j, eval( diff( R, j), params )],j in indets(R, name))];

Error, invalid input: eval expects its 2nd argument, eqns, to be of type {integer, equation, set(equation)}, but received [Lambda__p = 100, gamma__B = 0.5e-1, gamma__B = 0.5e-1, gamma__C = 0.1e-1, omega__C = 0.1e-2, omega__B = 0.1e-2, sigma__B = 0.1e-3, sigma__C = 0.1e-1, sigma__BC = 0.1e-1, theta__B = .8, theta__C = .5, mu__C = 1.0, Lambda__A = 1.0, Lambda__w = 1.0, 33*varphi__8, mu__A = 1.0, mu__w = 1.0, alpha__B = 0.5e-2, alpha__C = 0.5e-2, alpha__BC = .15, Zeta__B = .5, Zeta__C = .5, delta__A = .66, delta__w = 1.33]

 

#
# Compute all "sensitivities" (where the sensitivity
# is as defined in Rouben Rostamian response to the
# OP's earlier post) and evaluate numerically.
#
# For the purposes of this calculation "all" sensitivities
# means the sensitivity with respect to every variable
# returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( varName*diff( R, varName)/R, params )
# ]
#
  seq( [j, eval( j*diff( R, j)/R, params )],j in indets(R, name));

Error, invalid input: eval expects its 2nd argument, eqns, to be of type {integer, equation, set(equation)}, but received [Lambda__p = 100, gamma__B = 0.5e-1, gamma__B = 0.5e-1, gamma__C = 0.1e-1, omega__C = 0.1e-2, omega__B = 0.1e-2, sigma__B = 0.1e-3, sigma__C = 0.1e-1, sigma__BC = 0.1e-1, theta__B = .8, theta__C = .5, mu__C = 1.0, Lambda__A = 1.0, Lambda__w = 1.0, 33*varphi__8, mu__A = 1.0, mu__w = 1.0, alpha__B = 0.5e-2, alpha__C = 0.5e-2, alpha__BC = .15, Zeta__B = .5, Zeta__C = .5, delta__A = .66, delta__w = 1.33]

 

Download Computed_Sensitivity_Analys_for_CBD.mw

In the book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the another remarkable theorem was proved: any two flat bounded regions can be cut by a single straight line so that each of these regions is divided into two regions of equal area (the second  pancake problem). This is an existence theorem which does not provide any way to find this cut. In this post I made an attempt to find such cut for any 2 convex regions on the plane bounded by a piecewise smooth self-non-intersecting curves.
The Each_Into_2_Equal_Areas procedure returns a list of coordinates of 4 endpoints of the cutting segments. This procedure significantly uses my old procedures  Area  and  Picture , which can be found in detail at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  . The formal arguments of the Each_Into_2_Equal_Areas procedure are the lists  L1  and  L2 specifying the boundaries of the regions to be cut. When specifying  L1  and  L2 , the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing  t  with  -t  if necessary. The Pic procedure draws a picture of the source regions and cutting segments. For ease of use, the code for the  Area  and  Picture   procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

restart;
Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
abs(add(P[i], i = 1 .. nops(L))); 
end proc:

Picture := proc(L, C, N::posint := 100, Boundary::list := [linestyle = 1]) 
local i, var, var1, var2, e, e1, e2, P, Q, h; 
global Border;
uses plottools; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then P[i] := op(L[i]) else 
var := lhs(L[i, 2]); var1 := lhs(rhs(L[i, 2])); var2 := rhs(rhs(L[i, 2])); h := (var2-var1)/N; 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := seq(subs(var = var1+h*i, [e*cos(var), e*sin(var)]), i = 0 .. N) else 
P[i] := seq([var1+h*i, subs(var = var1+h*i, e)], i = 0 .. N) fi else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; P[i] := seq(subs(var = var1+h*i, [e1, e2]), i = 0 .. N) fi; fi; od; 
Q := [seq(P[i], i = 1 .. nops(L))]; Border := curve([op(Q), Q[1]], op(Boundary)); [polygon(Q, C), Border] 
end proc:

Each_Into_2_Equal_Areas:=proc(L1::list, L2::list)
local D, n, m, L10, L20, S1,S2, f, L11, L21, i, j, k, s, A, B, C , sol;

f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L10:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L1);
L20:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L2);
S1:=Area(L1); S2:=Area(L2);  
n:=nops(L1); m:=nops(L2);

for i from 1 to n do
for j from i to n do

for k from 1 to m do
for s from k to m do

if not ((nops({i,j})=1 and type(L1[i],listlist)) or (nops({k,s})=1 and type(L2[k],listlist))) then

A:=eval(L10[i,1],t=t1): 
B:=eval(L10[j,1],t=t2):
C:=eval(L20[k,1],t=t3): 
D:=eval(L20[s,1],t=t4):

L11:=`if`(j=i,[subsop([2,2]=t1..t2,L10[i]),[B,A]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]], [subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),op(L10[i+1..j-1]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]])):

L21:=`if`(s=k,[subsop([2,2]=t3..t4,L20[k]),[D,C]],`if`(s=k+1,[subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]], [subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),op(L20[k+1..s-1]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]])):

sol:=fsolve(simplify({Area(L11)-S1/2,Area(L21)-S2/2,eval(f(A,B),[x=C[1],y=C[2]]),eval(f(A,B),[x=D[1],y=D[2]])}),{t1=op([2,2,1],L10[i])..op([2,2,2],L10[i]),t2=op([2,2,1],L10[j])..op([2,2,2],L10[j]),t3=op([2,2,1],L20[k])..op([2,2,2],L20[k]),t4=op([2,2,1],L20[s])..op([2,2,2],L20[s])});

if type(sol,set(`=`)) then  return eval([A,B,C,D],sol) fi;

fi;
od: od: od: od:
end proc:

Pic := proc(L1,L2,col1,col2,Size:=[800,400])
local P1, P2, P3, T, P;
uses plots, plottools;
P1, P2 := Picture(L1, color=col1), Picture(L2, color=col2):
P3 := line(Sol[1..2][],color=red,thickness=3), line(Sol[3..4][],color=red,thickness=3), line(Sol[1],Sol[4],linestyle=2,thickness=3,color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
P:=pointplot(Sol, symbol=solidcircle, color=red, symbolsize=10);
display(P1,P2,P3,T,P, scaling=constrained, size=Size, axes=none);
end proc: 


Examples of use.

local D:
L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+16,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+16,5*sin(t)/2],t=Pi..2*Pi],[[21,0],[16,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue",[900,400]);

   

The specified regions may overlap:

L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+9,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+9,5*sin(t)/2],t=Pi..2*Pi],[[14,0],[9,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2):  Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

   


If there is a solution for which the cutting segments intersect the boundary of each of the regions at 2 points, then the procedure also works for such non-convex regions:

L1:=[[[cos(t),sin(t)],t=Pi/3..2*Pi-Pi/3],[[cos(-t)+1,sin(-t)],t=-Pi-Pi/3..-Pi+Pi/3]]:
L2:=[[[cos(t)+2,sin(t)],t=-Pi/6..Pi+Pi/6],[[cos(-t)+2,sin(-t)-1],t=-5*Pi/6..-Pi/6]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

   


A number of other interesting examples can be found in the attached file.

Each_Into_2_Equal_Areas1.mw

Hi,

I want to solve a first order ODE and (i) plot phi vs. x and then (ii) export data of the plotted curve as an T-shape ASCII file. I try to do this by Maple, but there are some errors and I couldn’t get the mentioned curve and ASCII file. Please, help me to remove errors:
(dphi/dx)**2+2*V(phi)=0,

where V(phi)= (1-alpha)*M^2 - (1-alpha)*M*sqrt(M^2 - 2*phi) + mu*(mu + beta*nu)*(1 - exp(phi/(mu + beta*nu))) + 
            (nu/beta)*(mu + beta*nu)*(1 - exp(beta*phi/(mu + beta*nu))) + (alpha/gamma1)*(1 - exp(-gamma1*phi)),

where we assume mu := 0.01, nu := 1 - mu, beta := 0.05, alpha := 0.3, M := sqrt(0.704), gamma1 := 0.001.
As it is seen from the attached figure, V(phi) and dV(phi)/dx =0 for phi=0 and phi=phi_m (two extreme points of V(phi)) and d^V(phi)/dphi^2<0 at phi=0 and phi=phi_m (phi and x are real). 

From the attached figure, it seems that phi=phi_m at x=0 and phi->0 as x-> infinity. (Hint: I expect the plot of V(phi) vs. phi and phi vs. x be similar to the curve "A" in the attached file).

Thanks.
rk4.mw

is(cos(x)=sqrt(1-sin(x)^2) or cos(x)=-sqrt(1-sin(x)^2)) assuming real;

                                          false  

In a simpler situation the result is correct

is(x=0 or x<>0) assuming real;

                                           true         

In the excellent book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the theorem is proved which states that any bounded planar region can be cut into 4 regions of equal area by 2 perpendicular cuts (the pancake problem). This is an existence theorem which does not provide any way to find these cuts. In this post I made an attempt to find such cuts for any convex region on the plane bounded by a piecewise smooth self-non-intersecting curve.
The Into_4_Equal_Areas procedure returns a list of coordinates of 5 points: the first 4 points are the endpoints of the cutting segments, the fifth point is the intersection point of these segments. This procedure significantly uses my old procedure Area , which can be found in detail at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  . The formal argument of the Into_4_Equal_Areas procedure is a list  L specifying the boundary of the region to be cut. When specifying  L, the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing  t  with  -t  if necessary. The Pic procedure draws a picture of the source region and cutting segments. For ease of use, the code for the  Area  procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

restart;
Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
abs(add(P[i], i = 1 .. nops(L))); 
end proc:

Into_4_Equal_Areas:=proc(L::list,N::symbol:='OneSolution', eps::numeric:=0)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, S, sol, Sol;
f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
S:=Area(L); c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
if not ((nops({i,j,k})=1 and type(L[i],listlist)) or (nops({j,k,m})=1 and type(L[j],listlist)))then
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
P:=eval([x,y], solve({f(A,C),f(B,D)},{x,y})):
L1:=`if`(j=i,[subsop([2,2]=t1..t2,L0[i]),[convert(B,list),P],[P,convert(A,list)]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]], [subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),L0[i+1..j-1][],subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]])):
L2:=`if`(k=j,[subsop([2,2]=t2..t3,L0[j]),[convert(C,list),P],[P,convert(B,list)]],`if`(k=j+1,[subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]], [subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),L0[j+1..k-1][],subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]])):
L3:=`if`(m=k,[subsop([2,2]=t3..t4,L0[k]),[convert(D,list),P],[P,convert(C,list)]],`if`(m=k+1,[subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]], [subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),L0[k+1..m-1][],subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]])):
sol:=fsolve({Area(L1)-S/4,Area(L2)-S/4,Area(L3)-S/4,LinearAlgebra:-DotProduct(D-B,C-A, conjugate=false)},{t1=op([2,2,1],L0[i])-eps..op([2,2,2],L0[i])+eps,t2=op([2,2,1],L0[j])-eps..op([2,2,2],L0[j])+eps,t3=op([2,2,1],L0[k])-eps..op([2,2,2],L0[k])+eps,t4=op([2,2,1],L0[m])-eps..op([2,2,2],L0[m])+eps}) assuming real:
if type(sol,set(`=`)) then if N='OneSolution' then return convert~(eval([A,B,C,D,P],sol),list) else c:=c+1; Sol[c]:=convert~(eval([A,B,C,D,P],sol),list) fi;
 fi; fi;
od: od: od: od:
convert(Sol,list);
end proc:

Pic:=proc(L,Sol)
local P1, P2, T;
uses plots, plottools;
P1:=seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L):
P2:=line(Sol[1],Sol[3],color=red, thickness=2), line(Sol[2],Sol[4],color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"],[Sol[5][],"P"]], font=[times,18], align=[left,above]);
display(P1,P2,T, scaling=constrained, size=[800,500], axes=none);
end proc: 


Examples of use:

L:=[[[0,0],[1,4]],[[1,4],[6,7]],[[6,7],[12,0]],[[12,0],[0,0]]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L, Sol);

# Check (areas of all 4 regions)
Area([[L[1,1],Sol[4],Sol[5],Sol[1],L[1,1]]]),
Area([[Sol[4],Sol[5],Sol[3],L[4,1],Sol[4]]]),
Area([[Sol[5],Sol[2],L[3,1],Sol[3],Sol[5]]]),
Area([[Sol[5],Sol[2],L[1,2],Sol[1],Sol[5]]]);

        


 

L:=[[[1+cos(-t),1+sin(-t)],t=-3*Pi/2..-Pi],[[0,1],[-1,0]],[[cos(t),sin(t)],t=Pi..2*Pi]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L,Sol);

    

# The boundary is the Archimedes spiral and the arc of a circle

L:=[[[t*cos(t),t*sin(t)],t=0..2*Pi],[[Pi+5*cos(-t),sqrt(25-Pi^2)+5*sin(-t)],t=arccos(Pi/5)..Pi-arccos(Pi/5)]]:
Sol:=evalf(Into_4_Equal_Areas(L));
Pic(L,Sol);

     

 

L:=[[[0,0],[2,0]],[[2,0],[1,sqrt(3)]],[[1,sqrt(3)],[0,0]]]:
Sol:=evalf[5](Into_4_Equal_Areas(L, AllSolutions, 0.1)); # All 3 solutions
plots:-display(<Pic(L, Sol[1]) | Pic(L, Sol[2])  | Pic(L, Sol[3])>, size=[300,300]);  


 

L:=[[[-t,-sin(-t)],t=-5*Pi/4..0],[[cos(t),sin(t)-1],t=Pi/2..3*Pi/2],[[t,cos(t)-3],t=0..3*Pi/2],[[3*Pi/2,-3],[5*Pi/4,sqrt(2)/2]]]:
Sol:=evalf(Into_4_Equal_Areas(L));
Pic(L,Sol);

More examples can be found in the attached file.

4_Equal_Area1.mw

[Edit]. The post has been edited. One inaccuracy in the code has been corrected, which could sometimes lead to errors. Two options have been added to the code of Into_4_Equal_Areas procedure. The first option is the formal argument  N . If N=OneSolution  (by default), the procedure returns one solution. If  N=AllSolutions , the procedure returns all solutions that it can find. The  eps  option has also been added (by default, eps=0). It is advisable to use it when we are looking for all solutions, and the ends of the cutting segments fall on the boundaries of intervals (this option slightly expands the boundaries of intervals, otherwise the  fsolve  command sometimes misses solutions). Two new examples have also been added.

 

Good afternoon.

I have a differential equation of non-integer degree and would like to know if it is possible to express a solution in terms of elementary or special-functions for certain values of the exponent, n>0.

For this equation, Maple provides an analytical solution for the exponent values n=0 and n=1, otherwise, there is no solution returned. I am particularly interested in the cases where n=1/2, 3/2, 2, 5/2, and 3

I am hoping that someone can help me resolve this - if a closed-form solution is not possible, then a numerical solution would also be welcome.

I have provided the details in the attached worksheet.

Thanks for reading!

MaplePrimes_Dec_19.mw

is(abs(x)=max(x,-x)) assuming real;

#  FAIL

I wonder if this will work in newer versions of Maple?

Hi,
How can I simplify this relation(See uploaded .mw file)?
For example, the second term is simplified as: 

deltae*(1-phi0/(kappa-3/2))^(-kappa+1/2)+(1/2)*deltab*(1-sqrt(2)*sqrt(1/(m*ub^2))*sqrt(-phi0));

di1.mw

How to insert a file inside another file, creating a session, without equations with the same name conflicting?

Hi,
I have an equation and I want to solve it parametrically to find x , but I couldn't do that with "solve" command. (I know x should be  real and positive). What should I do?
Root_of.mw

Hi. 

I am trying to solve a polynomial equation but the structure leads Maple to return a trivial solution and the other solutions are given as a RootOf expression. The equation involves a single variable, x, that is raised to a power, b and a multiplier, a (both are positive-valued). Please see attached worksheet.

I have not encountered this before and I cannot find a way to get to an explicit solution. Perhaps it is not possible (?).

Does anybody know how to deal with this? 

Thanks in advance ...

Roots_of_a_Polynomial_MaplePrimes.mw 

Good day, all.

I would like to explore the structure of the discrete modified form of the logistic equation.

In particular, I wish to plot the logistic-map to investigate the bifurcations of the system.

Is there a routine available in Maple that I can use?

I would like to consider the standard logistic equation with the inclusion of a shape parameter, m, introduced as a power law.

That is:

f(x) = a*x*(1-x^m)

where a > 0 denotes the growth rate, and m > 0  is a shape parameter. I wish to fix the value of a and take m to be the bifurcation parameter (so the logistic map would show m versus x for any given a).

Please note: The standard logistic equation (in discrete form) is given by f(x) = a*x*(1-x)

I would be grateful for any advice and support you can provide and I thank you for taking the time to read this.

Drear freinds,

I want to simplify f (a long experssion) in the form of f2. How to determine M1^2 and M0^2?

f1.mw

Hi,
How can I remove the mentioned error in attached worksheet?

s1.mw

1 2 3 4 5 6 7 Last Page 1 of 62