mmcdara

6468 Reputation

17 Badges

8 years, 43 days

MaplePrimes Activity


These are answers submitted by mmcdara

The first one consists in "rescaling the abscissas" of a classical ErrorPlot; the second one (which I prefer for it is more versatile) is to create "from scratch" your own error plot procedure.

ErrorPlot.mw

Example of what the second approach can produce:


As a rule I prefer  to define my own visualization procedure (of statistical data) instead of using those provided by the Statistics package.
This doesn't mean these latter are not "good" but only that they are not versatile enough to fit my requirements

Let us consider now this much simpler problem

(a*x+b*y)*diff(u(x, y), x) + (c*x+d*y)*diff(u(x, y), y) = 0

For this problem to have a solution you must have a=-d, which gives

                         1    2           1    2
               u(x, y) = - b y  + a x y - - c x 
                         2                2     

See compatibility.mw

So your problem could have a non trivial solution (aka a constant one) only for some values of values of alpha, b, c, d.
Here is an example of a non trivial solution for the case alpha=1, c=0, d=0 compatibility_2.mw

Another point of view                                                        
Your equation writes F(x, y, z, w) ° gradient(u(x, y, z, w) where F is field of directions in space [x, y, z, w] and ° denotes the dot product.
This direction field defines a family of lines such that the tangent to each of them is confounded with the direction of the field at the contact point.
Defining x(s), ..., w(s) as the cosine of this family gives:

                 d                              
                --- x(s) = alpha (b w(s) + y(s))
                 ds                             
                 d                              
                --- y(s) = -b w(s) + x(s) + z(s)
                 ds                             
                        d                
                       --- z(s) = -c y(s)
                        ds               
                    d                        
                   --- w(s) = d (y(s) - x(s))
                    ds                       

Integrating this system gives the expression of the characteristic lines.
Then u is an hypersurface made of these characteristic lines.

Look here  (characteristic system can be solved numerically but a numerical solution is provided as an illustration):

restart

E0 := alpha*(y+b*(w))*diff(u(x,y,z,w),x)+(x+z-b*(w))*diff(u(x,y,z,w),y)-c*y*diff(u(x,y,z,w),z)+d*(y-x)*diff(u(x,y,z,w),w)=0;

alpha*(b*w+y)*(diff(u(x, y, z, w), x))+(-b*w+x+z)*(diff(u(x, y, z, w), y))-c*y*(diff(u(x, y, z, w), z))+d*(y-x)*(diff(u(x, y, z, w), w)) = 0

(1)

eo := [op(lhs(E0))]

[alpha*(b*w+y)*(diff(u(x, y, z, w), x)), (-b*w+x+z)*(diff(u(x, y, z, w), y)), -c*y*(diff(u(x, y, z, w), z)), d*(y-x)*(diff(u(x, y, z, w), w))]

(2)

ed := map2(op, -1, eo)

[diff(u(x, y, z, w), x), diff(u(x, y, z, w), y), diff(u(x, y, z, w), z), diff(u(x, y, z, w), w)]

(3)

ec := eo /~ ed;

[alpha*(b*w+y), -b*w+x+z, -c*y, d*(y-x)]

(4)

ew := map2(op, -1, ed)

[x, y, z, w]

(5)

# Classical mathematical representation

zip((s1, s2) -> cat(d, s1)/s2 = ds, ew, ec)

[dx/(alpha*(b*w+y)) = ds, dy/(-b*w+x+z) = ds, -dz/(c*y) = ds, dw/(d*(y-x)) = ds]

(6)

# Maple translation

zip((s1, s2) -> diff(s1(s), s) = s2, ew, ec);

sys := lhs~(%) =~ eval(rhs~(%), ew =~ ew(s)):

print~(sys):

[diff(x(s), s) = alpha*(b*w+y), diff(y(s), s) = -b*w+x+z, diff(z(s), s) = -c*y, diff(w(s), s) = d*(y-x)]

 

diff(x(s), s) = alpha*(b*w(s)+y(s))

 

diff(y(s), s) = -b*w(s)+x(s)+z(s)

 

diff(z(s), s) = -c*y(s)

 

diff(w(s), s) = d*(y(s)-x(s))

(7)

Characteristic_lines := dsolve(sys, ew(s)):
whattype(%);
numelems(%%);

set

 

4

(8)

# Numerical appoximation

params := convert(indets(E0, name) minus convert(ew, set), list):

data := params =~ [seq(rand(0. .. 1.)(), i=1..numelems(params))];

[alpha = .6514625886, b = 0.2841433353e-1, c = 0.7936490070e-1, d = .3117315791]

(9)

numsys := eval({sys[], (ew(0) =~ [seq(rand(-1. .. 1.)(), i=1..4)])[]}, data)

{diff(w(s), s) = .3117315791*y(s)-.3117315791*x(s), diff(x(s), s) = 0.1851087527e-1*w(s)+.6514625886*y(s), diff(y(s), s) = -0.2841433353e-1*w(s)+x(s)+z(s), diff(z(s), s) = -0.7936490070e-1*y(s), w(0) = .442313422, x(0) = .359553901, y(0) = -.109372433, z(0) = 0.4759919e-2}

(10)

numsol := dsolve(numsys, 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) = 4, (2) = 4, (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.30045580019803245e-1, (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..4, {(1) = .442313422, (2) = .359553901, (3) = -.109372433, (4) = 0.4759919e-2}, 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..4, {(1) = .1, (2) = .1, (3) = .1, (4) = .1}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 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, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8]), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..4, {(1) = .442313422, (2) = .359553901, (3) = -.109372433, (4) = 0.4759919e-2}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = -.14617914657939401, (2) = -0.6306443973877118e-1, (3) = .35174577890249636, (4) = 0.8680332284362402e-2}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = w(s), Y[2] = x(s), Y[3] = y(s), Y[4] = z(s)]`; YP[1] := .3117315791*Y[3]-.3117315791*Y[2]; YP[2] := 0.1851087527e-1*Y[1]+.6514625886*Y[3]; YP[3] := -0.2841433353e-1*Y[1]+Y[2]+Y[4]; YP[4] := -0.7936490070e-1*Y[3]; 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] = w(s), Y[2] = x(s), Y[3] = y(s), Y[4] = z(s)]`; YP[1] := .3117315791*Y[3]-.3117315791*Y[2]; YP[2] := 0.1851087527e-1*Y[1]+.6514625886*Y[3]; YP[3] := -0.2841433353e-1*Y[1]+Y[2]+Y[4]; YP[4] := -0.7936490070e-1*Y[3]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..4, {(1) = 0., (2) = .442313422, (3) = .359553901, (4) = -.109372433}); _vmap := array( 1 .. 4, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4)  ] ); _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) = [s, w(s), x(s), y(s), z(s)], (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

(11)

with(plots):

col := [red, magenta, cyan, blue]:
display(
  seq(
    plots:-odeplot(numsol, [s, ew[i](s)], s=0..4, color=col[i], legend=ew[i](s))
    , i=1..4
  )
)

[red, magenta, cyan, blue]

 

 


A more general approach

AllParams := [params[],  S, seq(C__||i, i in ew) ];
numsys    := {sys[], seq(ew[i](S) = C__||(ew[i]), i=1..4)}

[alpha, b, c, d, S, C__x, C__y, C__z, C__w]

 

{diff(w(s), s) = d*(y(s)-x(s)), diff(x(s), s) = alpha*(b*w(s)+y(s)), diff(y(s), s) = -b*w(s)+x(s)+z(s), diff(z(s), s) = -c*y(s), w(S) = C__w, x(S) = C__x, y(S) = C__y, z(S) = C__z}

(12)

numsol := dsolve(numsys, numeric, parameters=AllParams)

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 := [alpha = alpha, b = b, c = c, d = d, S = S, C__x = C__x, C__y = C__y, C__z = C__z, C__w = C__w]; _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) = 4, (2) = 4, (3) = 0, (4) = 0, (5) = 9, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (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) = undefined, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = undefined, (6) = .0, (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..13, {(1) = C__w, (2) = C__x, (3) = C__y, (4) = C__z, (5) = Float(undefined), (6) = Float(undefined), (7) = Float(undefined), (8) = Float(undefined), (9) = Float(undefined), (10) = Float(undefined), (11) = Float(undefined), (12) = Float(undefined), (13) = Float(undefined)})), ( 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..4, {(1) = .1, (2) = .1, (3) = .1, (4) = .1}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 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, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8]), Array(1..13, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0}, datatype = float[8], order = C_order), Array(1..13, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0}, datatype = float[8], order = C_order), Array(1..13, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0}, datatype = float[8], order = C_order), Array(1..13, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..13, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0}, datatype = float[8], order = C_order), Array(1..13, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = w(s), Y[2] = x(s), Y[3] = y(s), Y[4] = z(s)]`; YP[1] := Y[8]*(Y[3]-Y[2]); YP[2] := Y[5]*(Y[1]*Y[6]+Y[3]); YP[3] := -Y[1]*Y[6]+Y[2]+Y[4]; YP[4] := -Y[7]*Y[3]; 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] = w(s), Y[2] = x(s), Y[3] = y(s), Y[4] = z(s)]`; YP[1] := Y[8]*(Y[3]-Y[2]); YP[2] := Y[5]*(Y[1]*Y[6]+Y[3]); YP[3] := -Y[1]*Y[6]+Y[2]+Y[4]; YP[4] := -Y[7]*Y[3]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..13, {(1) = S, (2) = C__w, (3) = C__x, (4) = C__y, (5) = C__z, (6) = undefined, (7) = undefined, (8) = undefined, (9) = undefined, (10) = undefined, (11) = undefined, (12) = undefined, (13) = undefined}); _vmap := array( 1 .. 4, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4)  ] ); _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) = [s, w(s), x(s), y(s), z(s)], (4) = [alpha = alpha, b = b, c = c, d = d, S = S, C__x = C__x, C__y = C__y, C__z = C__z, C__w = C__w]}); _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

(13)

# As a first illustration only the initial value S is changed while all the
# other parameters are ketp constant (same value 1, but you can change them)

V := proc(sigma)
  numsol(parameters=[1$4, sigma, 1$4]):
  display(
    seq(
      plots:-odeplot(numsol, [s, ew[i](s)], s=0..4, color=col[i], legend=ew[i](s))
      , i=1..4
    )
  )
end proc:

# Explore(V(sigma), parameters=[sigma=0. .. 10.])
 

# As a first illustration only the initial value S is changed while all the
# other parameters are ketp constant (same value 1, but you can change them)

V := proc(sigma)
  numsol(parameters=[1$4, sigma, 1$4]):
  odeplot(numsol, [s, x(s)], s=0..4, color=blue, legend=x(s))
end proc:

#display( seq(V(sigma), sigma=0..5, 0.1), view=[default, -2..2])

V := proc(sigma)
  numsol(parameters=[1$4, sigma, 1$4]):
  odeplot(numsol, [s, x(s), y(s)], s=0..4, color=blue)
end proc:

#display( seq(V(sigma), sigma=0..5, 0.1))

V := proc(sigma)
  numsol(parameters=[1$4, sigma, 1$4]):
  odeplot(numsol, [x(s), y(s), z(s)], s=0..4, color=blue)
end proc:

#display( seq(V(sigma), sigma=0..5, 0.1))

 

 

Download Characteristic_lines.mw

@Aixleft aimbot 

Here is a more readable code where some errors have been corrected.

Nevertheless it remains 4 expressions which "come from nowhere" and which I wasn't capable to check if the are correct or not: N11, N22, N33, N44 (it seems that you rewote yourself the expressions of N1, N2, N3, N4... Am I right?)

If it is so here are the expressions you are looking for:


                   2    1   /   (1/2)    \  (1/2)  3    / 4\
     tr :=    2 - v  + ---- \3 5      + 5/ 5      v  + O\v /
                  1200                                 

                       1   /   (1/2)    \  (1/2)  3    / 4\
     det :=       1 - ---- \7 5      + 5/ 5      v  + O\v /
                      1200                                 


What_I_understand.mw

@felixiao 

The output of taylor(...) is algebraic_expression+O(..) : the presence of O(..) causes the error you get.
The simplest thing to do is to use mtaylor instead whose output is an algebraic_expression.

A lengthy way is to write convert(taylor(...), polynom) instead of taylor(...).

MMSE_4COS_ok.mw

An advice for your future questions; avoid inserting the content of your uploaded file when it is that lond: insert only the link.

Illustrative example

restart:

PLSQ := (x, y) -> 1/x+1/y;

for i from -1 to 1 do
  for j from -1 to 1 do
    try
      print(i, j, PLSQ(i, j)):
    catch:
      print(i, j, "division by 0 not allowed")
    end try
  end do:
end do:
                           -1, -1, -2
               -1, 0, "division by 0 not allowed"
                            -1, 1, 0
               0, -1, "division by 0 not allowed"
               0, 0, "division by 0 not allowed"
               0, 1, "division by 0 not allowed"
                            1, -1, 0
               1, 0, "division by 0 not allowed"
                            1, 1, 2


to your last question  here
 

To know what are the methods Optimization uses, just do this

help("Methods Used by the Optimization Package")

More specifically, concerning Statistics:-NonlinearFit run this command

help("Regression Options")

From this help page:
 



The proof of the first formula is very simple but I rewrite F into F2 forreasons that will appear clearly during the the proof of your second formula.

Formula 1:

F := unapply(rsolve({F(1) = 1, F(2) = 1, F(n + 1) = F(n) + F(n - 1)}, F(n)), n):

rel2use := Phi = op([1, -1, 1], F(n)):
rel2use := { rel2use, C = op(1, F(n)) / rhs(rel2use)^n }

               /    1  (1/2)        1  (1/2)   1\ 
              { C = - 5     , Phi = - 5      + - }
               \    5               2          2/ 

# FORMULA 1

F2 := n -> C * (Phi^n -(1-Phi)^n):

ToProve_1 := F2(n + 1)^2 - F2(n)*F2(n+2) = (-1)^n:

expand~(ToProve_1):
simplify(eval(%, rel2use))
                             n       n
                         (-1)  = (-1) 


Formula 2:

# Preliminary:

G2 := n -> arccot(``(F2)(n));
              G2 := n -> arccot((F2)(n))

ToProve_2 := G2(2*n) - (G2(2*n + 1) + G2(2*n + 2)) = 0;
    arccot((F2)(2 n)) - arccot((F2)(2 n + 1)) - arccot((F2)(2 n + 2)) = 0

# From recurrence formula, replacing n+1 ny 2*n+2 and F by F2:

recur := ``(F2)(2*n+2) = ``(F2)(2*n+1) + ``(F2)(2*n);
           (F2)(2 n + 2) = (F2)(2 n + 1) + (F2)(2 n)

# One gets:

subs(recur, ToProve_2);
 arccot((F2)(2 n)) - arccot((F2)(2 n + 1)) - arccot((F2)(2 n + 1) + (F2)(2 n)) = 0

# Which has symbolic form:
# (U > 0, V > 0, U < V)

SymbolicForm := arccot(U) - arccot(V) - arccot(U+V)
             arccot(U) - arccot(V) - arccot(U + V)


# Proof

G2 := n -> arccot(F2(n)):

recur := F2(2*n+2) = F2(2*n+1) + F2(2*n):

ToProve_2 := simplify(G2(2*n) - (G2(2*n + 1) + G2(2*n + 2)) = 0):

SymbolicForm := -(1/2)*Pi + arctan(U) - arctan(V) - arctan(U+V)
           1                                           
         - - Pi + arctan(U) - arctan(V) - arctan(U + V)
           2                                           

# Thus 

SymbolicForm := -(1/2)*Pi + arctan((U-V)/(1+U*V)) - arctan(U+V):

# and next

SymbolicForm := simplify( -(1/2)*Pi + arctan(((U-V)/(1+U*V) - (U+V)) / (1+(U-V)/(1+U*V)*(U+V))) )
                              /  / 2          \ \
                 1            |V \U  + U V + 2/ |
               - - Pi - arctan|-----------------|
                 2            | 2          2    |
                              \U  + U V - V  + 1/

FocusOn := denom( op(1, indets(SymbolicForm, function)[]) )
                        2          2    
                       U  + U V - V  + 1

FocusOn := eval(FocusOn, {U=F2(2*n), V=F2(2*n+1)}):
FocusOn := eval(FocusOn, rel2use):
FocusOn := simplify(FocusOn) assuming n::integer;
                         (4 n)                     (4 n + 1)
      1  /  1  (1/2)   1\        1 /  1  (1/2)   1\         
    - -- |- - 5      + -|      + - |- - 5      + -|         
      10 \  2          2/        5 \  2          2/         

                            (4 n)       
         1  /  1  (1/2)   1\       (1/2)
       + -- |- - 5      + -|      5     
         10 \  2          2/            

FocusOn := ``(Y^(4*n)) * expand(eval(FocusOn, (-(1/2)*sqrt(5)+1/2)=Y) / Y^(4*n))
               / (4 n)\ /  1    1     1   (1/2)\
               \Y     / |- -- + - Y + -- 5     |
                        \  10   5     10       /

eval(FocusOn, Y=-(1/2)*sqrt(5)+1/2);
                               0
# Then 

'SymbolicForm' = op(1, SymbolicForm) + arctan(something/``(0))
                            1            /something\
           SymbolicForm = - - Pi + arctan|---------|
                            2            \   (0)   /

something := numer( op(1, indets(SymbolicForm, function)[]) )
                          / 2          \
                        V \U  + U V + 2/

# Obviously 'something' is strictly positive, thus 
#
# arctan(something/``(0)) = Pi/2;
#
# and thus

'ToProve_2' = ``(-Pi/2) + ``(Pi/2)

                           /  1   \    / 1   \     
               ToProve_2 = |- - Pi| + |  - Pi|
                           \  2   /    \ 2   /     

 

Download Fibonacci.mw


PS I believe the second proof could be easier to catch if one uses the trigonometric representation of the golden number.

@jalal 

" all calculators are programmed to divide by n... it's a broad debate"
Not all the calculators are programmed to divide by n (look to Excel(*) which proposes to divide by n or by n-1).
And it is not a broad debate: the situation where you have to divide by n is extremely clear, as is the situation where you have to divide by n-1.
The only situation when you have to divide by n is when your sample represents the whole population. If it is not the case dividing by n lead to a biased estimator of the variance, while dividing by n-1 lead to an unbiased estimator.

The reason I used n-1 (without asking you which one of these two situations you were in) is a direct consequence of the way you construct data and Weights:

  • As you randomly chose the elements of data from the range 4..20 this implicitely means that you draw a sample from a population whose indicuals are 4, 5, ..., 20. So data does not represents the whole population.
  • The same reasoning applies to Weights which are randomly chosen.

Had you said "I have population made of W[1] individuals X[1], ..., W[N] individuals X[N]" , the situation would have been different and the division by n perfectly justified.

 

 

Excel(*) : Excerpt from  Variance Excel:

               Excerpt from Population variance Excel

 

REMARK: Alternative ways to generate data are

N := 10:

data := sort(combinat:-randperm([$4..20])[1..N]);
data := sort(Statistics:-Shuffle([$2..20])[1..N]);


CodeTools:-Usage( sort(Statistics:-Shuffle([$2..20])[1..N]), iterations=10^5):
memory used=4.73KiB, alloc change=0 bytes, cpu time=46.96us, real time=47.03us, gc time=2.25us

CodeTools:-Usage( sort(combinat:-randperm([$4..20])[1..N]), iterations=10^5):
memory used=2.48KiB, alloc change=0 bytes, cpu time=15.79us, real time=15.80us, gc time=1.14us

 

@jalal 

...that you have weighted data and that you want to draw a boxplot of these data?

if it is so the trick is to build an unweighted sample from the data and their weights and draw the corresponding boxplot.
This is likely a hammer to kill a fly, but the method which consists in coding your own bowplot is far more complex (see Add-on 2 in attached file number two)


If the weights are strictly positive integers (which is the simpler case) use this:
 

restart

with(Statistics):

N       := 10:
data    := Sample(Normal(0, 1), N):

# non-weighted data

evalf[5](FivePointSummary(data));

Vector[column]([[minimum = -1.6349], [lowerhinge = -.61709], [median = 0.72465e-1], [upperhinge = 1.0065], [maximum = 1.7288]])

(1)

# weighted data
weights := [ seq(rand(1..10)(), n=1..N) ];

sample := `<,>`( seq(data[n]$weights[n], n=1..N) );
evalf[5](FivePointSummary(sample));

weights := [10, 2, 8, 10, 9, 1, 6, 7, 7, 3]

 

sample := Vector(4, {(1) = ` 1 .. 63 `*Vector[column], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

Vector(5, {(1) = minimum = -1.6349, (2) = lowerhinge = -.84476, (3) = median = -0.25428e-1, (4) = upperhinge = .21447, (5) = maximum = 1.7288})

(2)

BoxPlot([data, sample], axis[1]=[tickmarks=[1="unweighted", 2="weighted"]])

 

 


Download Boxplot_of_weighted_data.mw



If weights are strictly positive real numbers (a little bit more complex case) use this file
(look more specifically  the Maple issue of Add-on 3).

 

restart

with(Statistics):

Seed := 171803456472246;
randomize(Seed);

171803456472246

 

171803456472246

(1)

N       := 10:
data    := sort(Sample(Normal(0, 1), N)):
weights := [ seq(rand(0. .. 10. )(), n=1..N) ]:

# non-weighted data

evalf[5](FivePointSummary(data));
print():

interface(rtablesize=N+1):
VW := evalf[5](`<,>`(`<|>`("Values", "True weights"), `<|>`(data^+, < weights >)));

Vector[column]([[minimum = -2.3010], [lowerhinge = -.69574], [median = -.31568], [upperhinge = .55150], [maximum = 1.0886]])

 

NULL

 

VW := Matrix(11, 2, {(1, 1) = "Values", (1, 2) = "True weights", (2, 1) = -2.3010, (2, 2) = 4.3973, (3, 1) = -1.2398, (3, 2) = 9.3274, (4, 1) = -.69574, (4, 2) = 8.4820, (5, 1) = -.69352, (5, 2) = 1.0207, (6, 1) = -.52879, (6, 2) = 1.2482, (7, 1) = -.10258, (7, 2) = 2.8382, (8, 1) = .18670, (8, 2) = .62319, (9, 1) = .55150, (9, 2) = 5.5256, (10, 1) = .93235, (10, 2) = 5.5349, (11, 1) = 1.0886, (11, 2) = 8.4650})

(2)


FROM WEIGHTED TO UNWEIGHTED DATA

prob := weights /~ add(weights):
PT   := RandomVariable(ProbabilityTable(prob)):

UseHardwareFloats := false:

# The larger "LargeNumber", the closer the weights in sample to the real weights "weights"

LargeNumber := 10^5:
DataNumbers := round~(Sample(PT, LargeNumber)):

sample := Vector(LargeNumber, i -> data[DataNumbers[i]]):


evalf[5](FivePointSummary(sample));
print():

TS := Tally(sample):
ST := sort(evalf[4](lhs~(TS)) =~ rhs~(TS), key = (x -> lhs(x))):

# What are the weights we generate?

SW := evalf[5]( rhs~(ST) /~ (LargeNumber/add(weights)) ):

`<|>`(VW, Vector(["Weights in sample", SW[]]))
 

Vector[column]([[minimum = -2.3010], [lowerhinge = -1.2398], [median = -.52880], [upperhinge = .93235], [maximum = 1.0886]])

 

NULL

 

Matrix([["Values", "True weights", "Weights in sample"], [-2.3010, 4.3973, 4.5109], [-1.2398, 9.3274, 9.3545], [-.69574, 8.4820, 8.5130], [-.69352, 1.0207, .99388], [-.52879, 1.2482, 1.2245], [-.10258, 2.8382, 2.7960], [.18670, .62319, .62983], [.55150, 5.5256, 5.4635], [.93235, 5.5349, 5.4782], [1.0886, 8.4650, 8.4987]])

(3)

BoxPlot([data, sample], axis[1]=[tickmarks=[1="unweighted", 2="weighted"]])

 


Add-on 1

CONSTRUCTION OF AN UNWEIGHTED SAMPLE

IMO more complex than what I did above, but could be maybe an option when the size of data is large

Unweight := proc(Data, Weights, LargeNumber)
  local prob, PT, DataNumbers, sample:
   
  UseHardwareFloats := false:

  prob        := weights /~ add(weights):
  PT          := RandomVariable(ProbabilityTable(prob)):
  DataNumbers := round~(Sample(PT, LargeNumber)):
  sample      := Vector(LargeNumber, i -> data[DataNumbers[i]]):

  return sample
end proc:

BoxPlot([data, Unweight(data, weights, 10^5)], axis[1]=[tickmarks=[1="unweighted", 2="weighted"]])

 


Add-on 2

BOXPLOT RECONSTRUCTION OF WEIGHTED DATA

IMO more complex than what I did above, but could be maybe an option when the size of data is large

WBoxPlot := proc(Data, Weights, MyColor)
  local W_mean, W_median, W_Q25, W_Q75, W_deciles,
        W_upper, W_lower, W_outlier, W_boxplot,
        cumprob, k:

  uses plots, plottools:

  W_mean    := Mean    (Data, 'weights'=Weights):
  W_median  := Median  (Data, 'weights'=Weights):
  W_Q25     := Quantile(Data,  0.25, 'weights'=Weights):
  W_Q75     := Quantile(Data,  0.75, 'weights'=Weights):

  # This is the natural command you might thinkto to compute deciles.
  # But as shown in Add-on 3 it returns wrong results.
  # W_deciles := seq(Quantile(Data,  d, 'weights'=Weights), d=0.1..0.9, 0.1):
  #
  # So I use instead


  cumprob   := CumulativeSum(Weights) / add(Weights):
  W_deciles := NULL:
  for k from 0.1 to 0.9 by 0.1 do
    W_deciles := W_deciles, evalf[5]( data[ListTools:-BinaryPlace(cumprob, k)+1] )
  end do;

  W_upper   := min(max(Data), W_Q75+3/2*(W_Q75-W_Q25)):
  W_lower   := max(min(Data), W_Q25-3/2*(W_Q75-W_Q25)):
  W_outlier := Select((t -> is(t < Q_Lower or t > W_upper)), Data):

  W_boxplot := display(
    rectangle([-1, W_Q25], [1, W_Q75], color=MyColor)
    , plot([[-1, W_median], [1, W_median]], color=black)
    , plot([[0, W_Q75], [0, W_upper]], color=black)
    , plot([[0, W_Q25], [0, W_lower]], color=black)
    , plot([[-0.5, W_upper], [0.5, W_upper]], color=black)
    , plot([[-0.5, W_lower], [0.5, W_lower]], color=black)
    , pointplot([seq([0, W_deciles[i]], i=1..9)], symbol=circle, symbolsize=15, color=black)
  ):

  return W_boxplot
end proc:
 

use plots, plottools in
  display(
    BoxPlot([data, sample])
    , translate(scale(WBoxPlot(data, weights, "ForestGreen"), 0.75/2, 1), 3, 0)
    , axis[1]=[tickmarks=[1="unweighted", 2="weighted", 3="weighted like boxplot"]]
  )
end use

 


Add-on 3

A MAPLE ISSUE

Maple uses several methods to compute quantile (referenced by an integer from 1 to 8, see the corresponding help page).
The default method has number 7 and gives unexpected results when unequal weights are used.
This issue affects Percentile and Decile too.

So be extremely carefull while using these procedures with unequal weights

printf("This value is correct\n"):

Median    (data, 'weights'=weights);


printf("\nThese three values (method=default) are wrong\n"):

Quantile  (data,  0.5, 'weights'=weights);
Percentile(data,  50, 'weights'=weights);
Decile    (data,  5, 'weights'=weights);

# This is the reason why percentiles are not correctly placed on the green boxplot

This value is correct

 

-.5287852860

 


These three values (method=default) are wrong

 

-.2841943558

 

-.2841943558

 

-.2841943558

(4)

printf("\nThese three values (method=1) are correct\n"):

Quantile  (data,  0.5, 'weights'=weights, 'method'=1);
Percentile(data,  50, 'weights'=weights, 'method'=1);
Decile    (data,  5, 'weights'=weights, 'method'=1);
print():

printf("\nThese three values (method=2) are correct\n"):

Quantile  (data,  0.5, 'weights'=weights, 'method'=2);
Percentile(data,  50, 'weights'=weights, 'method'=2);
Decile    (data,  5, 'weights'=weights, 'method'=2);


These three values (method=1) are correct

 

-.5287852860

 

-.5287852860

 

-.5287852860

 

 


These three values (method=2) are correct

 

-.5287852860

 

-.5287852860

 

-.5287852860

(5)

# The computations are correct for the unweighted sample:

Median    (sample);
Quantile  (sample, 0.5);
Percentile(sample, 50);
Decile    (sample, 5);

-.5287852860

 

-.5287852860

 

-.5287852860

 

-.5287852860

(6)

# How to compute correctly deciles of weighted data

cumprob  := CumulativeSum(prob):

deciles := NULL:
for k from 0.1 to 0.9 by 0.1 do
  deciles := deciles, [k, evalf[5]( data[ListTools:-BinaryPlace(cumprob, k)+1] )]
end do:
`<,>`(< ':-Probability' | ':-Decile' >, convert([deciles], Matrix))

Matrix([[Probability, Decile], [.1, -1.2398], [.2, -1.2398], [.3, -.69574], [.4, -.69574], [.5, -.52879], [.6, .55150], [.7, .55150], [.8, .93235], [.9, 1.0886]])

(7)

# Check these values by using the unweighted sample

deciles := NULL:
for k from 1 to 9 do
  deciles := deciles, [k, evalf[5]( Decile(sample, k) )]
end do:
`<,>`(< ':-Probability' | ':-Decile' >, convert([deciles], Matrix))

Matrix([[Probability, Decile], [1, -1.2398], [2, -1.2398], [3, -.69574], [4, -.69574], [5, -.52879], [6, .55150], [7, .55150], [8, .93235], [9, 1.0886]])

(8)

# MAPLE BUG

wrong_deciles := NULL:
for k from 1 to 9 do
  wrong_deciles := wrong_deciles, [k, evalf[5]( Decile(data, k, 'weights'=prob) )]
end do:
`<,>`(< ':-Probability' | ':-Decile' >, convert([wrong_deciles], Matrix))

Matrix([[Probability, Decile], [1, -2.0633], [2, -1.5054], [3, -1.0750], [4, -.76037], [5, -.28440], [6, .38940], [7, .72001], [8, .96591], [9, 1.0564]])

(9)

# Check that BoxPlot correctlycomputes deciles for unweihghted sample

buw := BoxPlot(sample):
op(1..-2, select(has, [op(buw)], POINTS)[2]):
Deciles_from_BoxPlot := map(k -> [k, evalf[5](%[k][2])], [$1..9]):

`<,>`(< ':-Probability' | ':-Decile' >, convert(Deciles_from_BoxPlot, Matrix));

Matrix([[Probability, Decile], [1, -1.2398], [2, -1.2398], [3, -.69574], [4, -.69574], [5, -.52879], [6, .55150], [7, .55150], [8, .93235], [9, 1.0886]])

(10)

 

 

Download Boxplot_of_real_weighted_data.mw


NOTE:
There exist several quantile estimators for weighted data and the one I used here is just the most commonly used and easier to understand.

See for instance 
Akinshin


To moderators: please stop converting this comment into an answer, thank you


the expression you want to textplot is the result of the optimization of a Lenard-Jones model you use in your previous question?

if it is so the attached file contains a full reply (likely not the same data points the LJ model comes from).
 

restart;

pts := [[1.020408163, .9618194785], [1.021860259, .9591836735], [1.047329717, .9183673469], [1.077147904, .8775510204], [1.112217401, .8367346939], [1.153643677, .7959183673], [1.202786007, .7551020408], [1.224489796, .7394037725], [1.266731737, .7142857143], [1.350499968, .6734693878], [1.428571429, .6426647396], [1.463093584, .6326530612], [1.632653061, .5927592700], [1.638566951, .5918367347], [1.836734694, .5653094517], [2.024654644, .5510204082], [2.040816327, .5499100928], [2.244897959, .5415404104], [2.448979592, .5375577894], [2.653061224, .5364038532], [2.857142857, .5371174931], [3.061224490, .5390909888], [3.265306122, .5419306191], [3.469387755, .5453753581], [3.673469388, .5492483720], [3.759428515, .5510204082], [3.877551020, .5532080969], [4.081632653, .5571718264], [4.285714286, .5612391621], [4.489795918, .5653739340], [4.693877551, .5695506292], [4.897959184, .5737510616], [5.102040816, .5779621428], [5.306122449, .5821743709], [5.510204082, .5863807930], [5.714285714, .5905762856], [5.775983717, .5918367347], [5.918367347, .5943431266], [6.122448980, .5978993957], [6.326530612, .6014214453], [6.530612245, .6049104007], [6.734693878, .6083674439], [6.938775510, .6117937612], [7.142857143, .6151905094], [7.346938776, .6185587956], [7.551020408, .6218996660], [7.755102041, .6252141001], [7.959183673, .6285030098], [8.163265306, .6317672398], [8.219364040, .6326530612], [8.367346939, .6345736540], [8.571428571, .6371910423], [8.775510204, .6397830180], [8.979591837, .6423509843], [9.183673469, .6448962158], [9.387755102, .6474198713], [9.591836735, .6499230055], [9.795918367, .6524065794], [10.00000000, .6548714697]]:

model := a[1] + a[2]*((a[3]/Gamma)^a[4] - (a[5]/Gamma)^a[6]);

a[1]+a[2]*((a[3]/Gamma)^a[4]-(a[5]/Gamma)^a[6])

(1)

K := unapply(model, Gamma);
J := add((map2(op, 2, pts) -~ K~(map2(op, 1, pts)))^~2):
opt := evalf[6]( Optimization:-NLPSolve(J, assume=nonnegative) );

proc (Gamma) options operator, arrow; a[1]+a[2]*((a[3]/Gamma)^a[4]-(a[5]/Gamma)^a[6]) end proc

 

Warning, limiting number of major iterations has been reached

 

[0.810744e-4, [a[1] = .564441, a[2] = .350148, a[3] = 1.32558, a[4] = 2.94386, a[5] = 1.03995, a[6] = 1.69485]]

(2)

alias(STS = StringTools:-Substitute):


string_model := convert(model, string);
values       := evalf[3](opt[2])

"a[1]+a[2]*((a[3]/Gamma)^a[4]-(a[5]/Gamma)^a[6])"

 

[a[1] = .564, a[2] = .350, a[3] = 1.33, a[4] = 2.94, a[5] = 1.04, a[6] = 1.69]

(3)

for n from 1 to 6 do
  string_model := STS(string_model, convert(lhs(values[n]), string), convert(rhs(values[n]), string))
end do:

plots:-textplot([0, 0, cat("W__LJ = ", string_model)], axes=none,'font'=["helvetica","roman",25], size=[1000, 100])

 

 


 

Download display_formula_mmcdara.mw


(UPDATED: please do not turn this comment into an answer)


In your procedure the same quantity m is computed each time the loop is executed. So when theloop ends m contains the last result (here 11).
If you want to procedure to return all the values of m, you must record them in some structure.

For instance (tst returns a list)

tst := proc (M) 
  local i, m; 
  m := NULL; 
  for i from 2 to numelems(M)+1 do 
    m := m, M[1]-M[i-1] 
  end do; 
  return [m] 
end proc:


Or (tst returns a vector)

tst := proc (M) 
  local i, m; 
  m := Vector(numelems(M)):
  for i from 2 to numelems(M)+1 do 
    m[i-1] := M[1]-M[i-1] 
  end do; 
  return m
end proc:


Or more simply;

M[1] -~ M[2..-1]

 

@mz6687 

Here is the solution but it doesn't look to the image you sent before:
 

restart

with(LinearAlgebra):

A := kappa^3+(-(3*I)*m0-I*lambda-I*a[1]-I*a[2])*kappa^2+(-3*m0^2+(-2*lambda-2*a[1]-2*a[2])*m0+lambda^2-a[1]*a[2]+c^2)*kappa+I*m0^3+(I*lambda+I*a[1]+I*a[2])*m0^2+(-I*c^2-I*lambda^2+I*a[1]*a[2])*m0-I*lambda^3+(-I*a[1]-I*a[2])*lambda^2+(-I*c^2-I*a[1]*a[2])*lambda+(-I*c^2+I*c[2]^2)*a[2]-I*a[1]*c[2]^2:

B := (2*I)*c^6+(-(8*I)*c^2+(4*I)*a[1]^2+(4*I)*a[2]^2-(8*I)*c[1]^2-(8*I)*c[2]^2-(4*I)*n0+4*omega)*lambda^4+((8*I)*a[1]*c[1]^2+(8*I)*a[2]*c[2]^2)*lambda^3+(-(2*I)*c[2]^4+(-(6*I)*c^2+(2*I)*a[1]^2-(2*I)*c[1]^2-(4*I)*n0+4*omega)*c[2]^2+((2*I)*c^2-(2*I)*a[2]^2-(2*I)*c[1]^2)*a[1]^2+(2*I)*c^4+(-(8*I)*c[1]^2+(4*I)*n0-4*omega)*c^2+(4*I)*a[2]^2*c[1]^2-(2*I)*omega^2+(4*c[1]^2-4*n0)*omega+(2*I)*n0*(-2*c[1]^2+n0))*lambda^2+((4*I)*a[2]*c[2]^4+4*a[2]*(I*c^2-I*a[1]^2+I*c[1]^2+I*n0-omega)*c[2]^2+(8*(I*c^2-I*a[2]^2*(1/2)+I*n0*(1/2)-(1/2)*omega))*a[1]*c[1]^2)*lambda+(omega-I*c^2-I*a[2]^2-I*n0)*c[2]^4+((I*c^2+I*a[2]^2+I*c[1]^2+I*n0-omega)*a[1]^2-(2*I)*a[1]*a[2]*c[1]^2+(c^2-2*a[2]^2-c[1]^2)*(I*c^2+I*n0-omega))*c[2]^2+(-(2*I)*c^4+(I*a[2]^2-(2*I)*c[1]^2-(3*I)*n0+3*omega)*c^2+(I*c[1]^2+I*n0-omega)*a[2]^2+I*omega^2+omega*(c[1]^2+2*n0)-I*(c[1]^2+n0)*n0)*a[1]^2-(8*I)*lambda^6+(-I*a[2]^2+(5*I)*n0-5*omega)*c^4+((-(2*I)*n0+2*omega)*a[2]^2+(4*I)*n0^2-(4*I)*omega^2-8*n0*omega)*c^2+(-I*n0^2+I*omega^2+2*n0*omega)*a[2]^2-(3*I)*n0*omega^2+I*n0^3-3*n0^2*omega+omega^3:

x2 := simplify(coeff(A, kappa, 2)/I^(2-2+1)):

``

s3 := simplify(coeff(B, omega, 3)/I^(2-3+1)):

sol := solve({x2, x1, s2}, {m0, n0, lambda})

{lambda = -3*RootOf(12*_Z^2+(6*a[1]+6*a[2])*_Z+c^2+a[1]^2+a[1]*a[2]+a[2]^2)-a[1]-a[2], m0 = RootOf(12*_Z^2+(6*a[1]+6*a[2])*_Z+c^2+a[1]^2+a[1]*a[2]+a[2]^2), n0 = -(5/6)*c^2+(1/6)*a[1]^2-(5/6)*a[1]*a[2]-RootOf(12*_Z^2+(6*a[1]+6*a[2])*_Z+c^2+a[1]^2+a[1]*a[2]+a[2]^2)*a[1]+(1/6)*a[2]^2-RootOf(12*_Z^2+(6*a[1]+6*a[2])*_Z+c^2+a[1]^2+a[1]*a[2]+a[2]^2)*a[2]}

(1)

valsol := simplify~({allvalues(sol)});

{{lambda = -(1/4)*a[1]-(1/4)*a[2]-(1/4)*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2), m0 = -(1/4)*a[1]-(1/4)*a[2]+(1/12)*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2), n0 = -(5/6)*c^2+(5/12)*a[1]^2-(1/3)*a[1]*a[2]-(1/12)*a[1]*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)+(5/12)*a[2]^2-(1/12)*a[2]*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)}, {lambda = -(1/4)*a[1]-(1/4)*a[2]+(1/4)*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2), m0 = -(1/4)*a[1]-(1/4)*a[2]-(1/12)*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2), n0 = -(5/6)*c^2+(5/12)*a[1]^2-(1/3)*a[1]*a[2]+(1/12)*a[1]*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)+(5/12)*a[2]^2+(1/12)*a[2]*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)}}

(2)

nops(valsol);  # You gave two solutions, not one

2

(3)

# Here is for instance solution 1 (solution 2 has the same form)

print~(valsol[1]):

lambda = -(1/4)*a[1]-(1/4)*a[2]-(1/4)*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)

 

m0 = -(1/4)*a[1]-(1/4)*a[2]+(1/12)*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)

 

n0 = -(5/6)*c^2+(5/12)*a[1]^2-(1/3)*a[1]*a[2]-(1/12)*a[1]*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)+(5/12)*a[2]^2-(1/12)*a[2]*(-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)

(4)

# A simpler display

Z := 'Z':
R := select(type, indets(valsol), `^`)[] = Z:
ToDisplay := collect~(subs(R, valsol[1]), Z):

print~([(rhs=lhs)(R), ToDisplay[]]):

Z = (-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)

 

lambda = -(1/4)*a[1]-(1/4)*a[2]-(1/4)*Z

 

m0 = -(1/4)*a[1]-(1/4)*a[2]+(1/12)*Z

 

n0 = (-(1/12)*a[1]-(1/12)*a[2])*Z-(5/6)*c^2+(5/12)*a[1]^2-(1/3)*a[1]*a[2]+(5/12)*a[2]^2

(5)

ToDisplay

{lambda = -(1/4)*a[1]-(1/4)*a[2]-(1/4)*Z, m0 = -(1/4)*a[1]-(1/4)*a[2]+(1/12)*Z, n0 = (-(1/12)*a[1]-(1/12)*a[2])*Z-(5/6)*c^2+(5/12)*a[1]^2-(1/3)*a[1]*a[2]+(5/12)*a[2]^2}

(6)

# Or even simpler:

ToDisplay := map(u -> lhs(u) = value(``(numer(rhs(u))) / denom(rhs(u))), ToDisplay):
print~([(rhs=lhs)(R), ToDisplay[]]):

Z = (-12*c^2-3*a[1]^2+6*a[1]*a[2]-3*a[2]^2)^(1/2)

 

lambda = (1/4)*``(-a[1]-a[2]-Z)

 

m0 = (1/12)*``(-3*a[1]-3*a[2]+Z)

 

n0 = (1/12)*``(-Z*a[1]-Z*a[2]-10*c^2+5*a[1]^2-4*a[1]*a[2]+5*a[2]^2)

(7)

 

NULL


 

Download solA_mmcdara.mw


If you really want this "imaginary" form your imagedisplays (I don't find any interest in it), you can obtain it after some dymnastis:

eval(lhs(R), op(1, lhs(R)) = -op(1, lhs(R))*J^2):
newR := op(2, (combine(simplify(%), radical) assuming J > 0)) = Y:
eval( collect(subs(R, valsol[2][1]), Z), Z=J*Y):

s, r := selectremove(has, [op(rhs(%))], J):
convert(add(s) + J^2*(-add(r)), horner, J):
J * value(``(numer(op(1, %))) / denom(op(1, %))):

lambda = subs(J=I, %);
With = ( (rhs=lhs)(newR) )
                        1                        
               lambda = - I (I a[1] + I a[2] + Y)
                        4                        
         /                                             (1/2)\
         |    /    2         2                       2\     |
  With = \Y = \12 c  + 3 a[1]  - 6 a[1] a[2] + 3 a[2] /     /

 

@Andiguys 

# First thing: construct TRC correctly

TRC_ok := eval(TRC, {alpha = (u -> 0.02*exp(0.4*u)), beta = (u -> 0.2e-1*(1-exp(.5*I2)))})
                        12                   
          4.267000015 10   + 90 Q2 - 21.25 Q1

                             7               
             + 4.000000000 10  I1 exp(0.4 I1)

             + 2000000000 I2 (0.02 - 0.02 exp(0.5 I2))

# Assuming Q2 is nonnegative, TRC_ok is obviously mimimum when Q2=0.

TRC_reduced := eval(TRC_ok, Q2=0);
                12                            7               
  4.267000015 10   - 21.25 Q1 + 4.000000000 10  I1 exp(0.4 I1)

     + 2000000000 I2 (0.02 - 0.02 exp(0.5 I2))

# To be sure of what you are about to do verify what are the indeterminates of TRC_ok:

indets(TRC_ok, name)
                        {I1, I2, Q1, Q2}

# As you commented C1, C2 and C3, the Minimize command cannot work

C1 := 4e8 >= Q1,  Q1 >= .6e7;
C3 := I1 <= I2;
                              8      6      
                    Q1 <= 4 10 , 6 10  <= Q1
                            I1 <= I2

Next:

# Do you really want to write this:
#  C2 := if Q1= 4*(10)^(8) then Q2 =0 else Q2 = (x*d*h)-(delta*Q1*h)
# which meand that Q2 is always equal to (x*d-delta*Q1)*h excepted when Q2=4e8 where Q1=0 ???
#
# As you plane to use a numerical algorithm (Minimize) the probability that Q2=4e8 at some 
# iteration is exactly null... unless the minimum verifies Q1=4e-8 (limit of the first C1 
# condition).

Can we find "by hand" the minimum you are looking for?
TRC_reduced is obviously minimal when Q1 is maximal (Q1=4e8).

# The expression of TRC_reduced shows it is minimal when Q1 is maximal, thus 

evalf(eval(TRC_reduced, Q1=4e8));
               12         7                        9
       4.258 10   + 4.0 10  I1 exp(0.4 I1) + 2.0 10  I2 (0.02 - 0.02 exp(0.5 I2))

As a simple plot shows, the derivative wrt I2 of this expression is negative for 0 <= I2 <= 6.
Thus the minimum of TRC_reduced corresponds to Q1=4e8 and  I2=6.
It remains 

evalf(eval(TRC_reduced, [Q1=4e8, I2=6]));
                        12         7               
                4.254 10   + 4.0 10  I1 exp(0.4 I1)

... which is obviously minimum for I1=0 (a value which verifies I1 <= I2):

evalf(eval(TRC_reduced, [Q1=4e8, I2=6, I1=0]));
                             12
                     4.254 10  

(By the way I think you already asked  this same question and I already answered it. As I don't have time to verify this, please forgive me I'm wrong).

Whatever, here is a way to build an histogram from couples (bins, counts).
I assume your data are a list of equalities bin=count, where bin=a..b is a range and count is an integer.
For other structures, for instance a two-column matrix whose column 1 contains the bins and column 2 the counts, proceed as explained at the end oftheattached file.

Bins_and_counts_histogram.mw

I remain at your disposal for further help.

Remark: in the arrached file STEP 2 shows how to get the bins and counts of raw data while using TallyInto.

If you want to estimate the size and shape parameters of the Weibull random variable the data S could bea sample drawn from, you just have to do this

mle := MaximumLikelihoodEstimate(Weibull(a, b), S);
                 mle := [a = 129.944334168688, b = .9301095062]

If you want to build the corresponding random variable for any other purpose, I advice you to do this

W   := RandomVariable(Weibull(a, b)):
mle := MaximumLikelihoodEstimate(W, S):
W__opt := Specialize(W, mle):

# Example of use

N := numelems(S);
QuantilePlot(Sample(W__opt, N), S):
1 2 3 4 5 6 7 Last Page 1 of 55