Advanced Search

mmcdara

6149 Reputation

17 Badges

9 years, 62 days
Contact mmcdara

MaplePrimes Activity


These are answers submitted by mmcdara

Several errors fixed, the rest is up to ...

mmcdara 6149
June 02 2023
1 3


The solution is found, but the plots cannot be drawn because they still contain several quantities that have not been evaluated:
 

``

> 

restart; with(plots)
> 

PDEtools[declare]((F, T, G, H)(Y), prime = Y):

F(Y)*`will now be displayed as`*F

  

T(Y)*`will now be displayed as`*T

  

G(Y)*`will now be displayed as`*G

  

H(Y)*`will now be displayed as`*H

  

`derivatives with respect to`*Y*`of functions of one variable will now be displayed with '`

 (1)
> 

p1 := 0.1e-1:
> 

rf := 1050:
> 

sigma1 := 25000:
> 

sigma2 := 0.210e-5:
> 

sigma3 := 6.30*10^7:
> 

sigma4 := 10^(-10):
> 

sigma5 := 1.69*10^7:
> 

sigma6 := 4.10*10^7:
> 

NULL
> 

M := 1:
> 

alp := .1:
> 

NULL
> 

NULL
> 

B1 := 1+2.5*p+6.2*p^2:
> 

a1 := B1*p1+B2*p2+B3*p3:
> 

a2 := 1-p1-p2-p3+p1*rs1/rf+p2*rs2/rf+p3*rs3/rf:
> 

a3 := 1-p1-p2-p3+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf)+p3*rs3*cps3/(rf*cpf):
> 

a4 := B4*p1+B5*p2+B6*p3:
> 

NULL
> 

a5 := 1+3*((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3)/(2+(p1*sigma1+p2*sigma2+p3*sigma3)/((p1+p2+p3)*sigmaf)-((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3)):
> 

NULL
> 

> 

a6 := B1*p1+B2*p2+B3*p3:
> 

a7 := 1-p1-p2-p3+p1*rs4/rf+p2*rs5/rf+p3*rs6/rf:
> 

a8 := 1-p1-p2-p3+p1*rs4*cps4/(rf*cpf)+p2*rs5*cps5/(rf*cpf)+p3*rs6*cps6/(rf*cpf):
> 

a9 := B7*p1+B8*p2+B9*p3:
> 

NULL
> 

a10 := 1+3*((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3)/(2+(p1*sigma4+p2*sigma5+p3*sigma6)/((p1+p2+p3)*sigmaf)-((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3)):
> 

W := sum(b[i]*Y^i, i = 0 .. 3);

Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0]

  

Y^3*c[3]+Y^2*c[2]+Y*c[1]+c[0]

  

Y^2*d[2]+Y*d[1]+d[0]

  

Y^2*h[2]+Y*h[1]+h[0]

 (2)
> 

F := a1*(1+1/bet)*(diff(W, `$`(Y, 2)))+a2*Ra*(diff(W, Y))+A-a5*M*W-S2*W^2+a2*Gr*Theta-S*betu*(W-U) = 0;

9.1682928*Y*b[3]+3.0560976*b[2]+2.433571428*Y^2*b[3]+1.622380952*Y*b[2]+.8111904760*b[1]+1-1.346703274*b[3]*Y^3-1.346703274*b[2]*Y^2-1.346703274*b[1]*Y-1.346703274*b[0]-.5*(Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+.8111904760*c[3]*Y^3+.8111904760*c[2]*Y^2+.8111904760*c[1]*Y+.8111904760*c[0]+0.5e-1*d[2]*Y^2+0.5e-1*d[1]*Y+0.5e-1*d[0] = 0

 (3)
> 

T := (a4+Rd)*(diff(Theta, `$`(Y, 2)))+a3*Pr*Ra*(diff(Theta, Y))+Q*Theta+Pr*alp*S*bett*(Theta-Phi)+Pr*Ec*((1+1/bet)*a1*(diff(W, Y))^2+a5*M*W^2+(1+1/bet)*a1*S1*W^2+S2*W^3+S*betu*(W-U)) = 0;

-.525*d[0]-.105*h[0]+.205*c[0]+10.23967791*c[1]+2.266560888*c[2]+.525*b[0]+6.799682664*Y*c[3]+30.71903373*Y^2*c[3]+20.47935582*Y*c[2]-.105*Y^2*h[2]-.105*Y*h[1]+.525*b[1]*Y+.525*b[2]*Y^2+.525*b[3]*Y^3+.205*c[1]*Y+.205*c[2]*Y^2+.205*c[3]*Y^3-.525*d[1]*Y-.525*d[2]*Y^2+21.63764058*(Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+16.04451240*(3*Y^2*b[3]+2*Y*b[2]+b[1])^2+5.25*(Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^3 = 0

 (4)
> 

G := Ra*(diff(U, Y))+betu*(W-U) = 0;

1.0*Y*d[2]+.5*d[1]+.5*b[3]*Y^3+.5*b[2]*Y^2-.5*d[2]*Y^2+.5*b[1]*Y-.5*d[1]*Y+.5*b[0]-.5*d[0] = 0

 (5)
> 

H := Ra*(diff(Phi, Y))+bett*(Theta-Phi) = 0;

1.0*Y*h[2]+.5*h[1]+.5*c[3]*Y^3+.5*c[2]*Y^2-.5*Y^2*h[2]+.5*c[1]*Y-.5*Y*h[1]+.5*c[0]-.5*h[0] = 0

 (6)
> 

BCS := (D(W))(0) = 0, (D(Theta))(0) = 0, W(1) = -delta*(1+1/bet)*(D(W))(1), Theta(1) = 1+g*(D(Theta))(1), U(1) = -delta*(1+1/bet)*(D(W))(1), Phi(1) = 1+g*(D(Theta))(1):
> 

> 

W     := unapply(W, Y);
F     := unapply(F, Y);
U     := unapply(U, Y);
G     := unapply(G, Y);
Phi   := unapply(Phi, Y);
Theta := unapply(Theta, Y);
T     := unapply(T, Y);
H     := unapply(H, Y);

proc (Y) options operator, arrow; Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0] end proc

  

proc (Y) options operator, arrow; 9.1682928*Y*b[3]+3.0560976*b[2]+2.433571428*Y^2*b[3]+1.622380952*Y*b[2]+.8111904760*b[1]+1-1.346703274*Y^3*b[3]-1.346703274*Y^2*b[2]-1.346703274*Y*b[1]-1.346703274*b[0]-.5*(Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+.8111904760*c[3]*Y^3+.8111904760*c[2]*Y^2+.8111904760*c[1]*Y+.8111904760*c[0]+0.5e-1*d[2]*Y^2+0.5e-1*d[1]*Y+0.5e-1*d[0] = 0 end proc

  

proc (Y) options operator, arrow; d[2]*Y^2+d[1]*Y+d[0] end proc

  

proc (Y) options operator, arrow; 1.0*Y*d[2]+.5*d[1]+.5*Y^3*b[3]+.5*Y^2*b[2]-.5*d[2]*Y^2+.5*Y*b[1]-.5*d[1]*Y+.5*b[0]-.5*d[0] = 0 end proc

  

proc (Y) options operator, arrow; Y^2*h[2]+Y*h[1]+h[0] end proc

  

proc (Y) options operator, arrow; c[3]*Y^3+c[2]*Y^2+c[1]*Y+c[0] end proc

  

proc (Y) options operator, arrow; -.525*d[0]-.105*h[0]+.205*c[0]+10.23967791*c[1]+2.266560888*c[2]+.525*b[0]+6.799682664*Y*c[3]+30.71903373*Y^2*c[3]+20.47935582*Y*c[2]-.105*Y^2*h[2]-.105*Y*h[1]+.525*Y*b[1]+.525*Y^2*b[2]+.525*Y^3*b[3]+.205*c[1]*Y+.205*c[2]*Y^2+.205*c[3]*Y^3-.525*d[1]*Y-.525*d[2]*Y^2+21.63764058*(Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+16.04451240*(3*Y^2*b[3]+2*Y*b[2]+b[1])^2+5.25*(Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^3 = 0 end proc

  

proc (Y) options operator, arrow; 1.0*Y*h[2]+.5*h[1]+.5*c[3]*Y^3+.5*c[2]*Y^2-.5*Y^2*h[2]+.5*c[1]*Y-.5*Y*h[1]+.5*c[0]-.5*h[0] = 0 end proc

 (7)
> 

 

> 

z1 := (D(W))(0) = 0:
> 

z2 := (D(Theta))(0) = 0:
> 

z3 := W(1) = -delta*(1+1/bet)*(D(W))(1):
> 

z4 := Theta(1) = 1+g*(D(Theta))(1):
> 

z5 := U(1) = -delta*(1+1/bet)*(D(W))(1):
> 

z6 := Phi(1) = 1+g*(D(Theta))(1):
> 

z7 := F(0):
> 

z8 := T(0):
> 

z9 := G(0):
> 

z10 := H(0):
> 

z11 := F(1):
> 

z12 := T(1):
> 

z13 := G(1):
> 

z14 := H(1):
> 

Zs :=  [z1, z2, z3, z4, z5, z6, z7, z8, z9, z10, z11, z12, z13, z14]:
print~(Zs):

b[1] = 0

  

c[1] = 0

  

b[3]+b[2]+b[1]+b[0] = -0.6e-2*b[3]-0.4e-2*b[2]-0.2e-2*b[1]

  

c[3]+c[2]+c[1]+c[0] = 1+.3*c[3]+.2*c[2]+.1*c[1]

  

d[2]+d[1]+d[0] = -0.6e-2*b[3]-0.4e-2*b[2]-0.2e-2*b[1]

  

h[2]+h[1]+h[0] = 1+.3*c[3]+.2*c[2]+.1*c[1]

  

1.+3.0560976*b[2]+.8111904760*b[1]-1.346703274*b[0]-.5*b[0]^2+.8111904760*c[0]+0.5e-1*d[0] = 0

  

-.525*d[0]-.105*h[0]+.205*c[0]+10.23967791*c[1]+2.266560888*c[2]+.525*b[0]+21.63764058*b[0]^2+16.04451240*b[1]^2+5.25*b[0]^3 = 0

  

.5*d[1]+.5*b[0]-.5*d[0] = 0

  

.5*h[1]+.5*c[0]-.5*h[0] = 0

  

10.25516095*b[3]+3.331775278*b[2]-.5355127980*b[1]+1-1.346703274*b[0]-.5*(b[3]+b[2]+b[1]+b[0])^2+.8111904760*c[3]+.8111904760*c[2]+.8111904760*c[1]+.8111904760*c[0]+0.5e-1*d[2]+0.5e-1*d[1]+0.5e-1*d[0] = 0

  

-.525*d[0]-.105*h[0]+.205*c[0]+10.44467791*c[1]+22.95091671*c[2]+.525*b[0]+37.72371639*c[3]-.105*h[2]-.105*h[1]+.525*b[1]+.525*b[2]+.525*b[3]-.525*d[1]-.525*d[2]+21.63764058*(b[3]+b[2]+b[1]+b[0])^2+16.04451240*(3*b[3]+2*b[2]+b[1])^2+5.25*(b[3]+b[2]+b[1]+b[0])^3 = 0

  

.5*d[2]+.5*b[3]+.5*b[2]+.5*b[1]+.5*b[0]-.5*d[0] = 0

  

.5*h[2]+.5*c[3]+.5*c[2]+.5*c[1]+.5*c[0]-.5*h[0] = 0

 (8)
> 

IZs := indets(Zs);

numelems(Zs), numelems(IZs);

{b[0], b[1], b[2], b[3], c[0], c[1], c[2], c[3], d[0], d[1], d[2], h[0], h[1], h[2]}

  

14, 14

 (9)
> 

Zsol := fsolve(Zs)

{b[0] = .6178879176, b[1] = 0., b[2] = -.7968855995, b[3] = .1810986325, c[0] = 3.017782504, c[1] = 0., c[2] = -4.488567511, c[3] = 2.247245007, d[0] = .2073632729, d[1] = -.4105246447, d[2] = .2052623223, h[0] = 1.523567501, h[1] = -1.494215002, h[2] = .7471075012}

 (10)
> 

# Note that some b, c, d, h are still unknown
# Which means plots will not be obtained

F := unapply(eval(sum(b[i]*Y^i,i=0..5), Zsol), Y);
T := unapply(eval(sum(c[i]*Y^i,i=0..5), Zsol), Y);
G := unapply(eval(sum(d[i]*Y^i,i=0..5), Zsol), Y);
H := unapply(eval(sum(h[i]*Y^i,i=0..5), Zsol), Y);

proc (Y) options operator, arrow; b[5]*Y^5+b[4]*Y^4+.1810986325*Y^3-.7968855995*Y^2+.6178879176 end proc

  

proc (Y) options operator, arrow; c[5]*Y^5+c[4]*Y^4+2.247245007*Y^3-4.488567511*Y^2+3.017782504 end proc

  

proc (Y) options operator, arrow; d[5]*Y^5+d[4]*Y^4+d[3]*Y^3+.2052623223*Y^2-.4105246447*Y+.2073632729 end proc

  

proc (Y) options operator, arrow; h[5]*Y^5+h[4]*Y^4+h[3]*Y^3+.7471075012*Y^2-1.494215002*Y+1.523567501 end proc

 (11)
> 

``
> 

plot(F(Y), Y = 0 .. 1, axes = boxed, color = green, thickness = 2, labels = [Y, W]);

Warning, expecting only range variable Y in expression b[5]*Y^5+b[4]*Y^4+.1810986325*Y^3-.7968855995*Y^2+.6178879176 to be plotted but found names [b[4], b[5]]

  

 
> 

# Maybe the range i=0..5 should be replaced by i=0..3 ?

F := unapply(eval(sum(b[i]*Y^i,i=0..3), Zsol), Y);
plot(F(Y), Y = 0 .. 1, axes = boxed, color = green, thickness = 2, labels = [Y, W])

proc (Y) options operator, arrow; .1810986325*Y^3-.7968855995*Y^2+.6178879176 end proc

  

 
> 

``


 

Download AGM_mmcdara.mw

Unnecessary complexity...

mmcdara 6149
June 01 2023
0 1

and a lot of typos (x and y instead of X and Y for instance).
More of this, a formal solution does exist

> 

restart:
> 

with(plots):
> 

eq1 := diff(X(t), t) - Y(t) = 0;
eq2 := diff(Y(t), t) + X(t) = 0;
bcs := X(0) = 1, Y(0) = 0;

diff(X(t), t)-Y(t) = 0

  

diff(Y(t), t)+X(t) = 0

  

X(0) = 1, Y(0) = 0

 (1)
> 

sol := dsolve({bcs, eq1, eq2})

{X(t) = cos(t), Y(t) = -sin(t)}

 (2)
> 

plot(eval([X(t), Y(t)], sol), t=0..10, color=[green, blue], legend=[X(t), Y(t)])

 
> 

solnum := dsolve({bcs, eq1, eq2}, numeric);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.5047658755841546e-2, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..2, {(1) = 1.0, (2) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = 1.0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = -1.0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = X(t), Y[2] = Y(t)]`; YP[2] := -Y[1]; YP[1] := Y[2]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = X(t), Y[2] = Y(t)]`; YP[2] := -Y[1]; YP[1] := Y[2]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..2, {(1) = 0., (2) = 1.}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, X(t), Y(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (3)
> 

plots:-odeplot(solnum, [[t, X(t)], [t, Y(t)]], t=0..10, color=[green, blue], legend=[X(t), Y(t)])

 
> 

 

Download Excat_and_Num_Sol.mw

REMARK: using declare doesn't authorize you to forget the independent variable t:

> 

restart
> 

with(PDETools, declare):
> 

declare(X(t)); declare(Y(t))
> 

eq1 := diff(X(t), t)-Y(t) = 0;

diff(X(t), t)-Y(t) = 0

  

diff(Y(t), t)+X(t) = 0

  

X(0) = 1, Y(0) = 0

  

{X(t) = cos(t), Y(t) = -sin(t)}

 (1)
> 

``

Download Declare_exact_sol.mw

@ogunmiloro An initial condition is...

mmcdara 6149
May 31 2023
1 0

@ogunmiloro 

An initial condition is not just a value but this value PLUS  the time at wihich this value holds.
I'm waiting for values for a and b in the equality

G[dp](a) = b


You say F[d]=1.00... is this "1" or "1 million"  (same thing for L[d] and E[c])?

Last point: as I told you before let C the defined by the relation 

C := ( psi[1]*1.00 + psi[2]*2. 018 +psi[3]*1.10 ) *10^6

Any triple (psi[1],  psi[2], psi[3]) which keeps C constant gives exactly the same solution G[dp](T): thus psi[1],  psi[2], psi[3] are not identifiable and your problem depends only on two parameters: C and delta[3].


Taking arbitrarily the initial condition

G[dp](2015) = 2.11634848e6;

One get these results
 

> 

restart


PARTIAL DATA (year 2023 removed for reasons explained below)

> 

data := Matrix(
  8, 2,
  [
    2022, 3023744.14,
    2021, 2804897.12,
    2020, 2707438.88,
    2019, 2969316.28,
    2018, 2734527.11,
    2017, 2590856.01,
    2016, 2445951.47,
    2015, 2116348.48
  ]
);


FLE    := [F[d](T) = 1e6, L[d](T) = 2.018103e6, E[c](T) = 1.1e6];
edo    := diff(G[dp](T), T) = C-delta[3]*G[dp](T);
ic     := G[dp](2015) = 2.11634848e6;
sys    := eval({edo, ic}, FLE)

data := Matrix(8, 2, {(1, 1) = 2022, (1, 2) = 3023744.14, (2, 1) = 2021, (2, 2) = 2804897.12, (3, 1) = 2020, (3, 2) = 2707438.88, (4, 1) = 2019, (4, 2) = 2969316.28, (5, 1) = 2018, (5, 2) = 2734527.11, (6, 1) = 2017, (6, 2) = 2590856.01, (7, 1) = 2016, (7, 2) = 2445951.47, (8, 1) = 2015, (8, 2) = 2116348.48})

  

[F[d](T) = 0.1e7, L[d](T) = 2018103., E[c](T) = 0.11e7]

  

diff(G[dp](T), T) = C-delta[3]*G[dp](T)

  

G[dp](2015) = 2116348.48

  

{diff(G[dp](T), T) = C-delta[3]*G[dp](T), G[dp](2015) = 2116348.48}

 (1)
> 

sol := rhs(dsolve(sys))

C/delta[3]-exp(-delta[3]*T)*(-52908712/25+C/delta[3])/exp(-2015*delta[3])

 (2)


Define the objective function as the (square of the) ℓ2 norm of data[.., 2] minus the predictions

Firstly set c = C/delta[3] to avoid singulatities.

> 

eval(sol, C =~ c*~delta[3]):
simplify(%) assuming delta[3] > 0:
SOL := unapply(%, T);

pred := simplify~(SOL~(data[..,1])):

ObjFunc := evalf(add((data[..,2]-pred)^~2));

proc (T) options operator, arrow; c+(52908712/25)*exp(-delta[3]*(T-2015))-exp(-delta[3]*(T-2015))*c end proc

  

(3023744.14-1.*c-2116348.480*exp(-7.*delta[3])+exp(-7.*delta[3])*c)^2+(2804897.12-1.*c-2116348.480*exp(-6.*delta[3])+exp(-6.*delta[3])*c)^2+(2707438.88-1.*c-2116348.480*exp(-5.*delta[3])+exp(-5.*delta[3])*c)^2+(2969316.28-1.*c-2116348.480*exp(-4.*delta[3])+exp(-4.*delta[3])*c)^2+(2734527.11-1.*c-2116348.480*exp(-3.*delta[3])+exp(-3.*delta[3])*c)^2+(2590856.01-1.*c-2116348.480*exp(-2.*delta[3])+exp(-2.*delta[3])*c)^2+(2445951.47-1.*c-2116348.480*exp(-1.*delta[3])+exp(-1.*delta[3])*c)^2

 (3)


With Optimization:-Minimize

> 

Opt := Optimization:-Minimize(ObjFunc, initialpoint={delta[3] = 1, c=1e6}, assume=nonnegative);

[69400773501.3077698, [c = HFloat(2938688.7509472966), delta[3] = HFloat(0.47654470252163217)]]

 (4)
> 

OptPred := eval(pred, Opt[2]):
plots:-display(
  dataplot(data[.., 1], data[.., 2], color=blue, legend="data"),
  dataplot(data[.., 1], OptPred, color=red, legend="Best pred.")
)

 


By solving grad(ObjFunc) = 0

> 

eqs := {diff(ObjFunc, c), diff(ObjFunc, delta[3])}:
> 

h := 5;
extremum := fsolve(eqs, {c=1e6, delta[3]=1});

5

  

{c = 2938688.896, delta[3] = .4765445243}

 (5)
> 

# check if this extremum is a minimum (partial derivatives must be < 0)

eval(ObjFunc, extremum);
eval(eqs, extremum);

0.6940077346e11

  

{-0.199e4, -0.293e-2}

 (6)
> 

BestPred := evalf(eval(pred, extremum)):

plots:-display(
  dataplot(data[.., 1], data[.., 2], color=blue, legend="data"),
  dataplot(data[.., 1], BestPred, color=red, legend="Best pred.")
)

 


COMPLETE DATA (year 2023 included)

Note :the value of the GDP at year 2023 being extremely small, it is likely that the model won't fit the data in a correct way

> 

data := Matrix(
  9, 2,
  [
    2023, 1764580,
    2022, 3023744.14,
    2021, 2804897.12,
    2020, 2707438.88,
    2019, 2969316.28,
    2018, 2734527.11,
    2017, 2590856.01,
    2016, 2445951.47,
    2015, 2116348.48
  ]
):


Define the objective function as the (square of the) ℓ2 norm of data[.., 2] minus the predictions

Firstly set c = C/delta[3] to avoid singulatities.

> 

pred := simplify~(SOL~(data[..,1])):

ObjFunc := evalf(add((data[..,2]-pred)^~2)):


With Optimization:-Minimize

> 

Opt := Optimization:-Minimize(ObjFunc, initialpoint={delta[3] = 1, c=1e6}, assume=nonnegative);

[1066972164798.57336, [c = HFloat(2661106.488854825), delta[3] = HFloat(1.197406240485557)]]

 (7)
> 

OptPred := eval(pred, Opt[2]):
plots:-display(
  dataplot(data[.., 1], data[.., 2], color=blue, legend="data"),
  dataplot(data[.., 1], OptPred, color=red, legend="Best pred.")
)

 
> 

 

 

Download fitting_mmcdara.mw

The value of  C is c*delta[3] (= 3.186425515*10^6 for the last result), but there are an infinity of triples (psi[1],  psi[2], psi[3]) which give this same value.


                                              ADD-ON
In case F[d], L[d] and E[c] are unknown functions of the time whose only values at (let's say) the end of each year are known you can proceed this way (two implicit discretization scheme of your ode are avaliable).
Other approaches can be used, such as filtering for instance.
 

> 

restart


PARTIAL DATA (year 2023 removed , see previous file)

Rows are reversed for future use.

> 

data := Matrix(
  8, 2,
  [
    2022, 3023744.14,
    2021, 2804897.12,
    2020, 2707438.88,
    2019, 2969316.28,
    2018, 2734527.11,
    2017, 2590856.01,
    2016, 2445951.47,
    2015, 2116348.48
  ]
):
data := data[[seq](8..1, -1)]

data := Matrix(8, 2, {(1, 1) = 2015, (1, 2) = 2116348.48, (2, 1) = 2016, (2, 2) = 2445951.47, (3, 1) = 2017, (3, 2) = 2590856.01, (4, 1) = 2018, (4, 2) = 2734527.11, (5, 1) = 2019, (5, 2) = 2969316.28, (6, 1) = 2020, (6, 2) = 2707438.88, (7, 1) = 2021, (7, 2) = 2804897.12, (8, 1) = 2022, (8, 2) = 3023744.14})

 (1)
> 

# Simulation of varions of Fd, Ld and Ec from year to year.
# For the sake of simplicity I used random variations here.

sigma := 1e4;
Fd := Statistics:-Sample(Normal(1e6, sigma), 8)^+:
Ld := Statistics:-Sample(Normal(2.018103e6, sigma), 8)^+:
Ec := Statistics:-Sample(Normal(1.1e6, sigma), 8)^+:

0.1e5

 (2)
> 

# Examples of numerical schemes (time step h = 1 year)

(g(x+h)-g(x))/h = (psi__1*F(x)+psi__2*L(x)+psi__3*E(x))-delta__3*(g(x)+g(x+h))/2;
isolate(%, g(x+h)):
expand(eval(%, h=1));


# With lumping

(g(x+h)-g(x))/h = (psi__1*(F(x)+F(x+h))/2+psi__2*(L(x)+L(x+h))/2+psi__3*(E(x)+E(x+h))/2)-delta__3*(g(x)+g(x+h))/2;
isolate(%, g(x+h)):
expand(eval(%, h=1));

(g(x+h)-g(x))/h = psi__1*F(x)+psi__2*L(x)+psi__3*E(x)-(1/2)*delta__3*(g(x)+g(x+h))

  

g(x+1) = 2*psi__1*F(x)/(delta__3+2)+2*psi__2*L(x)/(delta__3+2)+2*psi__3*E(x)/(delta__3+2)-delta__3*g(x)/(delta__3+2)+2*g(x)/(delta__3+2)

  

(g(x+h)-g(x))/h = (1/2)*psi__1*(F(x)+F(x+h))+(1/2)*psi__2*(L(x)+L(x+h))+(1/2)*psi__3*(E(x)+E(x+h))-(1/2)*delta__3*(g(x)+g(x+h))

  

g(x+1) = psi__1*F(x)/(delta__3+2)+psi__1*F(x+1)/(delta__3+2)+psi__2*L(x)/(delta__3+2)+psi__2*L(x+1)/(delta__3+2)+psi__3*E(x)/(delta__3+2)+psi__3*E(x+1)/(delta__3+2)-delta__3*g(x)/(delta__3+2)+2*g(x)/(delta__3+2)

 (3)
> 

# Use of a

Predictor := proc(p1, p2, p3, d)
  local N, Gdp, i:
  N      := numelems(data[..,1]);
  Gdp    := Vector(N):
  Gdp[1] := rhs(ic):
  for i from 2 to N do
  # Gdp[i] := (Gdp[i] + p1*Fd[i] + p2*Ld[i] + p3*Ec[i] )/(1+d/2) - Gdp[i-1]*d/2/(1+d/2):
    Gdp[i] := (Gdp[i] + p1*(Fd[i]+Fd[i-1])/2 + p2*(Ld[i]+Ld[i-1])/2 + p3*(Ec[i]+Ec[i-1])/2)/(1+d/2) - Gdp[i-1]*d/2/(1+d/2):
  end do:
  return Gdp              
end proc:

# Objective function to minimize (RSS=Residual Sum of Squares).

RSS := proc(p1, p2, p3, d)
  local N, Gdp, P, i:
  N      := numelems(data[..,1]);
  Gdp    := Vector(N):
  Gdp[1] := rhs(ic):
  for i from 2 to N do
  # Gdp[i] := (Gdp[i] + p1*Fd[i] + p2*Ld[i] + p3*Ec[i])/(1+d/2) - Gdp[i-1]*d/2/(1+d/2):
    Gdp[i] := (Gdp[i] + p1*(Fd[i]+Fd[i-1])/2 + p2*(Ld[i]+Ld[i-1])/2 + p3*(Ec[i]+Ec[i-1])/2)/(1+d/2) - Gdp[i-1]*d/2/(1+d/2):
  end do:

  evalf(add((data[..,2]-~Gdp)^~2));
end proc:
> 

# Initial condition (notional example).

ic := G[dp](2015) = 2.11634848e6;

G[dp](2015) = 2116348.48

 (4)


With Optimization:-Minimize

Two different initial points used.

 

> 

omega := 0:
Opt1 := Optimization:-Minimize(RSS, initialpoint=[0.7$3, 0.1]);
Opt2 := Optimization:-Minimize(RSS, initialpoint=[0.1$4]);

Opt1 := [4.95342259811396484*10^11, Vector(4, {(1) = -18.09280538479193, (2) = 80.83226827021075, (3) = -119.56599085212119, (4) = 4.864058833001402}, datatype = float[8])]

  

Opt2 := [1.24601467797302109*10^11, Vector(4, {(1) = 3.3486953190379873, (2) = .10770503061396366, (3) = -2.9747682860150406, (4) = -.46363134805073025}, datatype = float[8])]

 (5)
> 

# Synthesis & plots

OptPred1 := Predictor(entries(Opt1[2], nolist)):
OptPred2 := Predictor(entries(Opt2[2], nolist)):
`<,>`(
  `<|>`(["Year", "GDP", "Prediction 1", "Prediction 2"]),
  evalf[3](`<|>`(data, OptPred1, OptPred2))
);

plots:-display(
  dataplot(data[.., 1], data[.., 2], color=blue, legend="data"),
  dataplot(data[.., 1], OptPred1, color=red, legend="Best pred. 1"),
  dataplot(data[.., 1], OptPred2, color=black, legend="Best pred. 2")
)

Matrix([["Year", "GDP", "Prediction 1", "Prediction 2"], [0.202e4, 0.212e7, 0.212e7, 0.212e7], [0.202e4, 0.245e7, 0.286e7, 0.231e7], [0.202e4, 0.259e7, 0.271e7, 0.275e7], [0.202e4, 0.273e7, 0.250e7, 0.285e7], [0.202e4, 0.297e7, 0.265e7, 0.283e7], [0.202e4, 0.271e7, 0.287e7, 0.275e7], [0.202e4, 0.280e7, 0.274e7, 0.276e7], [0.202e4, 0.302e7, 0.267e7, 0.281e7]])

  

 
> 

 

 


Download fitting_mmcdara_2.mw

If you want to inforce positivity of the fitting parameters do this

Opt1 := Optimization:-Minimize(RSS, initialpoint=[0.1$3, 0.1], assume=nonnegative);
OptPred1 := Predictor(entries(Opt1[2], nolist)):
`<,>`(
  `<|>`(["Year", "GDP", "Prediction 1"]),
  evalf[3](`<|>`(data, OptPred1))
);

plots:-display(
  dataplot(data[.., 1], data[.., 2], color=blue, legend="data"),
  dataplot(data[.., 1], OptPred1, color=red, legend="Best pred. 1")
)

Watch out: remain extremely careful about the solution you get and check that some components are not stuck on the initial point.

Could you check if...

mmcdara 6149
May 30 2023
2 0

if arrow procedure Nb_characters works correctly (the test strings have been generate with a French Mac keyboard)

Download Maybe_this.mw

@AHSAN Have you take the time to lo...

mmcdara 6149
May 28 2023
0 5

@AHSAN 

Have you take the time to look here https://www.mapleprimes.com/posts/210159-2D-Contour-Plot-And-Legend ?
I think it would be a good starting point

If I'm not mistaken...

mmcdara 6149
May 28 2023
0 0


Based upon Descarte's tule of signs,  and assuming gamma > 0, I would say there is only one positive root.
If gamma < 0 I can' be that conclusive.
Only_1_positive_root.mw

A solution...

mmcdara 6149
May 28 2023
1 3

Here is a solution (an example about how to place the legend within the plot and to build a composite legend made of line style and symbol ifs proposed in the file):
A_solution.mw

Download A_solution.mw

You can also refer to the answer I gave to your earlier question question https://www.mapleprimes.com/questions/236446-How-To-To-Draw-The-Horizontal-And-Veritcal

Here are two ways:...

mmcdara 6149
May 27 2023
1 6

The first option is immediate, the second one places t "legend-ike box" within the plot window itself.

> 

restart:
> 

N := (-18*sqrt(2)*(M-4/3)*(x^2+2)^2*arctan((1/2)*x*sqrt(2))-9*Pi*(M-4/3)*(x^2+2)^2*sqrt(2)+(-36*M+48)*x^3+(-120*M+96)*x)/(4*(x^2+2)^2)

(1/4)*(-18*2^(1/2)*(M-4/3)*(x^2+2)^2*arctan((1/2)*x*2^(1/2))-9*Pi*(M-4/3)*(x^2+2)^2*2^(1/2)+(-36*M+48)*x^3+(-120*M+96)*x)/(x^2+2)^2

 (1)
> 

p := plot(
  [subs(M = 1.3015, N)]
  , x = -5 .. 5
  , labels = ["x", "N"]
  , axes = boxed
  , labeldirections = ["horizontal", "horizontal"]
  , colour = red
  , linestyle = solid
  , legend = ["N= 1, M= 1.3015"]
  , legendstyle = [location = bottom, font = ["TIMES", "italic", 15]]
 ):

ymin, ymax := (min, max)(plottools:-getdata(p)[3][.., 2]):
eps   := 0.1:
Vline := plot([[0, ymin-eps], [0, ymax+eps]], color=black):
Hline := plot([[-5, 0], [5, 0]], color=black):

plots:-display(p, Vline, Hline);

 
> 

p := plot(
  [subs(M = 1.3015, N)]
  , x = -5 .. 5
  , labels = ["x", "N"]
  , axes = boxed
  , labeldirections = ["horizontal", "horizontal"]
  , colour = red
  , linestyle = solid
  , gridlines = false
 ):


ymin, ymax := (min, max)(plottools:-getdata(p)[3][.., 2]):
eps   := 0.1:
Vline := plot([[0, ymin-eps], [0, ymax+eps]], color=black):
Hline := plot([[-5, 0], [5, 0]], color=black):

Left  := -4.8:
Right := -0.5:

Legend := plots:-display(
  plot(ymin, x=Left..Left+1, color=red)
  , plots:-textplot([Left+1, ymin, " N= 1, M= 1.3015"], align=right)
  , plottools:-rectangle([-4.9, ymin-eps], [-0.2, ymin+eps], color=white)
  #, plot([[Left-0.1, ymin-eps], [Right, ymin-eps], [Right, ymin+eps], [Left-0.1, ymin+eps], [Left-0.1, ymin-eps]], color=black)
):

plots:-display(p, Vline, Hline, Legend);

 
> 

 

Download TwoOptions.mw

Waiting for concurrent proposals......

mmcdara 6149
May 27 2023
0 1

... here is a first answer.
The coding is simple and probably not that efficient vor (very) large matrices, Let's just say it's the start of the exchanges to come?

> 

restart
> 

with(LinearAlgebra):
> 

M := RandomMatrix(6$2, generator=0..10)

M := Matrix(6, 6, {(1, 1) = 5, (1, 2) = 6, (1, 3) = 6, (1, 4) = 1, (1, 5) = 1, (1, 6) = 5, (2, 1) = 10, (2, 2) = 10, (2, 3) = 5, (2, 4) = 0, (2, 5) = 2, (2, 6) = 4, (3, 1) = 1, (3, 2) = 5, (3, 3) = 3, (3, 4) = 3, (3, 5) = 4, (3, 6) = 7, (4, 1) = 2, (4, 2) = 1, (4, 3) = 6, (4, 4) = 9, (4, 5) = 3, (4, 6) = 8, (5, 1) = 0, (5, 2) = 3, (5, 3) = 5, (5, 4) = 8, (5, 5) = 6, (5, 6) = 3, (6, 1) = 6, (6, 2) = 4, (6, 3) = 0, (6, 4) = 3, (6, 5) = 9, (6, 6) = 4})

 (1)
> 

# A rough code: it"s likely that adjustments must be done for large matrices
# in order to get an efficient code.

SmallestBelow := proc(M, withdiag)
  local coding, minval, here:
  coding := x -> piecewise(x = 0, 1, 0);
  map(x -> [lhs(x)], op(coding~(M))[3]);
  if withdiag then
    minval := min(seq(seq(M[i, j], j=1..i), i=1..numelems(M[1])));
    here   := select((x -> is(x[2] <= x[2])), map(x -> [lhs(x)], op(coding~(M))[3]));
  else
    minval := min(seq(seq(M[i, j], j=2..i), i=1..numelems(M[1])));
    here   := select((x -> is(x[2] < x[1])), map(x -> [lhs(x)], op(coding~(M))[3]));
  end if:
  minval, here
end proc:

# Strictly delow the diagonal

minval, here := SmallestBelow(M, false);


# Diagonal included

minval, here := SmallestBelow(M, true);

0, {[5, 1], [6, 3]}

  

0, {[2, 4], [5, 1], [6, 3]}

 (2)
> 

 

Download Min_value_and_locations.mw

Obvious solution:...

mmcdara 6149
May 26 2023
0 0

(x-4*y)^9+3*z reaches its maximum value when z=17 and when (x-4*y)^9 is maximum, thus when x=3 and y=0.
Maximum value = 3^9+3*17 = 19734

Didn't I already answered  simi...

mmcdara 6149
May 25 2023
1 3

Didn't I already answered a similar question weeks ago ???

In case you didn't read my answer, which can be possible, here it is again.
In case you read it but didn't like it, I would have been happy to receive your comments instead of seeing it again under hidden form.
By the way, the code below works for any version of Maple.

> 

restart;
> 

#with(student):
> 

U := a[0] + sum(-a[i]*tanh(xi[n])^i, i = 1 .. 1) + sum(-b[i]*tanh(xi[n])^(-i), i = 1 .. 1):
> 

u(xi[n + 1]) := a[0] - a[1]*(tanh(xi[n]) + tanh(d))/(1 + tanh(xi[n])*tanh(d)) - b[1]*(1 + tanh(xi[n])*tanh(d))/(tanh(xi[n]) + tanh(d)):
> 

u(xi[n - 1]) := a[0] - a[1]*(tanh(xi[n]) - tanh(d))/(1 - tanh(xi[n])*tanh(d)) - b[1]*(1 - tanh(xi[n])*tanh(d))/(tanh(xi[n]) - tanh(d)):
> 

eq := c[1]^2*diff(U, xi[n], xi[n])*ln(1 + U/alpha) - beta*(u(xi[n + 1]) - 2*U - u(xi[n - 1])):
> 

fin1 := simplify(numer(eq)):
> 

fin := simplify(subs(tanh(xi[n]) = Psi, fin1)):
> 

# As the degree(Psi*ln(Psi), Psi) returns FAIL, the idea is to replace ln(...) by
# something else, let's say H.
#
# Observe that EQ[1] is  the identity 0=0

lnPsi := select(has, indets(fin, function), ln)[]:
fin_H := eval(fin, lnPsi=H):
deg   := degree(fin_H, Psi);
print():

for i from 0 to deg do
  collect(coeff(fin_H, Psi, i), H);
  EQ[i] := eval(%, H=lnPsi) = 0;
end do:

print~(EQ):

10

  

 

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*b[1]*c[1]^2 = 0

  

0 = 0

  

(2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2+2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*b[1] = 0

  

2*tanh(d)^3*beta*a[1]+2*tanh(d)^2*beta*a[0]+2*tanh(d)*beta*b[1] = 0

  

(-2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2-2*tanh(d)^2*b[1]*c[1]^2-2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*b[1]-2*tanh(d)^2*beta*a[1]+2*beta*b[1] = 0

  

-2*beta*a[0]-2*tanh(d)^4*beta*a[0]-2*tanh(d)^3*beta*a[1]-2*tanh(d)^3*beta*b[1]-2*tanh(d)*beta*a[1]-2*tanh(d)*beta*b[1] = 0

  

(-2*tanh(d)^4*a[1]*c[1]^2-2*tanh(d)^2*a[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2-2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*a[1]-2*tanh(d)^2*beta*b[1]+2*beta*a[1] = 0

  

2*tanh(d)^2*beta*a[0]+2*tanh(d)*beta*a[1]+2*tanh(d)^3*beta*b[1] = 0

  

0 = 0

  

(2*tanh(d)^4*a[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2+2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*a[1] = 0

  

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*a[1]*c[1]^2 = 0

 (1)
> 

 

Download Degree_3_mmcdara.mw

@Kitonum's answer is likely the simp...

mmcdara 6149
May 23 2023
0 0

@Kitonum's answer is likely the simpler and the most efficient.

Nevertheless here is a slight variant:

for i from 3 to nops(n) do 
  if i::odd then
    f[i] := n[i-2]-3*n[i-1]+3*n[i]; 
  else
    f[i] := n[i-3]+4*n[i-2]-6*n[i-1]+n[i]; 
  end if:
end do:

@Rouben Rostamian   Hi Rouben, Not...

mmcdara 6149
May 22 2023
2 2

@Rouben Rostamian  

Hi Rouben,
Not in my intention to interfere in your exchange with the OP, but I suggest you to do something like this, whith your more recent version than my old 2015.2 (I was forced to rotate the plot to circumvent the lack of the orientation option in texplot)

 

> 

restart;
> 

f := theta -> 6000*sin(theta) - 3819.718635

proc (theta) options operator, arrow; 6000*sin(theta)-3819.718635 end proc

 (1)
> 

pl1 := plots:-display(
    plot(f(t)+3819.718635, t=0..Pi, coords=polar, color=red, thickness=3)
  , seq(plot(r, t=0..Pi, coords=polar, color=gray, thickness=1), r in [seq](0..7000, 3819.718635/5))
  , plot(3819.718635, t=0..Pi, coords=polar, color=black, thickness=1)
  , seq(plot([[0, 0], [7000*cos(t*Pi), 7000*sin(t*Pi)]], color=gray), t in [seq](0..1, 1/6))
  , axes=none
  , scaling=constrained
):

# In Maple the option orientation is doesn't exist to oriente the text.
# So this trick to suggest you a modification of your answer.

plots:-display(
    plottools:-rotate(pl1, Pi/2)
   , seq(plots:-textplot([0, r, nprintf("%1.2e", r-3819.718635)], align=right), r in [seq](0..6500, 3819.718635/5))
   , seq(
       plots:-textplot(
         [
           7000*cos(t*Pi+Pi/2)
           , 7000*sin(t*Pi+Pi/2)
           , typeset(Pi-t*Pi)
         ]
         , align={left, piecewise(t=1/2, NULL, t > 1/2, below, above)}
       )
       , t in [seq](0..1, 1/6)
     )
);

 
> 

 


 

Download suggestion.mw

@delvin Is it this that you want?&n...

mmcdara 6149
May 21 2023
0 1

@delvin 

Is it this that you want?

> 

restart;
> 

#with(student):
> 

U := a[0] + sum(-a[i]*tanh(xi[n])^i, i = 1 .. 1) + sum(-b[i]*tanh(xi[n])^(-i), i = 1 .. 1):
> 

u(xi[n + 1]) := a[0] - a[1]*(tanh(xi[n]) + tanh(d))/(1 + tanh(xi[n])*tanh(d)) - b[1]*(1 + tanh(xi[n])*tanh(d))/(tanh(xi[n]) + tanh(d)):
> 

u(xi[n - 1]) := a[0] - a[1]*(tanh(xi[n]) - tanh(d))/(1 - tanh(xi[n])*tanh(d)) - b[1]*(1 - tanh(xi[n])*tanh(d))/(tanh(xi[n]) - tanh(d)):
> 

eq := c[1]^2*diff(U, xi[n], xi[n])*ln(1 + U/alpha) - beta*(u(xi[n + 1]) - 2*U - u(xi[n - 1])):
> 

fin1 := simplify(numer(eq)):
> 

fin := simplify(subs(tanh(xi[n]) = Psi, fin1)):
> 

# As degree(Psi*ln(Psi), Psi) returns FAIL, the idea is to replace ln(...) by
# something else, let's say H.
#
# Observe that EQ[1] is  the identity 0=0

lnPsi := select(has, indets(fin, function), ln)[]:
fin_H := eval(fin, lnPsi=H):
deg   := degree(fin_H, Psi);
print():

for i from 0 to deg do
  collect(coeff(fin_H, Psi, i), H);
  EQ[i] := eval(%, H=lnPsi) = 0;
end do:

print~(EQ):

10

  

 

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*b[1]*c[1]^2 = 0

  

0 = 0

  

(2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2+2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*b[1] = 0

  

2*tanh(d)^3*beta*a[1]+2*tanh(d)^2*beta*a[0]+2*tanh(d)*beta*b[1] = 0

  

(-2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2-2*tanh(d)^2*b[1]*c[1]^2-2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*b[1]-2*tanh(d)^2*beta*a[1]+2*beta*b[1] = 0

  

-2*beta*a[0]-2*tanh(d)^4*beta*a[0]-2*tanh(d)^3*beta*a[1]-2*tanh(d)^3*beta*b[1]-2*tanh(d)*beta*a[1]-2*tanh(d)*beta*b[1] = 0

  

(-2*tanh(d)^4*a[1]*c[1]^2-2*tanh(d)^2*a[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2-2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*a[1]-2*tanh(d)^2*beta*b[1]+2*beta*a[1] = 0

  

2*tanh(d)^2*beta*a[0]+2*tanh(d)*beta*a[1]+2*tanh(d)^3*beta*b[1] = 0

  

0 = 0

  

(2*tanh(d)^4*a[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2+2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*a[1] = 0

  

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*a[1]*c[1]^2 = 0

 (1)
> 

 

Download Degree_3_mmcdara.mw

Here is a simpler and corrected version ...

mmcdara 6149
May 16 2023
0 6

Here is a simpler and corrected version of your code: note the Alert Maplet displays the content of the TextArea component and no longer "test".
But the three times knocking is still here
SimplerButStill3Times.mw

I will keep digging to try and understand where it comes from but I'm afraid it could be some undesired "side effect" between DocumentTools and Maplets packages (?).

This suspicion comes from this observation:
Here is a slight variant of the previous file
SimplerButStill3Times_2.mw
Proceed this way:

  1. write a number in the text area
  2. click on the shortcut to execute the "line" you wrote
  3. click on the line "void"

The  Alert Maplet is displayed once
Changing the number and redoing steps 2 and 3 still displays the Alert Maplet once

First 8 9 10 11 12 13 14 Last Page 10 of 52