mmcdara

Ask a Question

Create a Post

18 Badges

Member for: 8 years, 129 days

Contact mmcdara

MaplePrimes Activity

These are answers submitted by mmcdara

I hope this will enlighten you enough to...

Answered: mmcdara 6635

February 20 2024

1 0

@jalal

To get this histogram you present:

Given the nature of your data I suppose they are already consist in a serie of elements of the type [C[k], N[k]] where
- C[k] is the range of the kth class,
- N[k] is the number of elements C[k] contains.

This is what I call synthetic data by opposition to raw data (see below).

If it is so you have two possibilities:
1. To build an Statistics:-Histogram you must have raw data, not synthetic data.
  So either you have these raw data or you don't.
  - If you have them, let R the list/vector of incomes, just use Statistics:-Histogram with the ad hoc customization (first plot in the attached file).
  - If you do not have them, you can recreate fake raw data from synthetic data so that their histogram is the one you want (second plot in the attached file).
2. If you only have synthetic data and bo not necessarily want to use Statistics:-Histogram, then you can proceed as explaine in the third plot in the attached file.

restart

CASE 1: YOU HAVE RAW DATA

Simulation of hypothetic raw data

>	R := Statistics:-Sample(Triangular(6, 97, 53), 10^3):

Raw data histogram

C := [seq(0..100, 10)]:
Statistics:-Histogram(
  R
  , frequencyscale=absolute
  , binbounds=C
  , axis[1]=[tickmarks=C]
  , view=[-1..101, default]
  , color="Pink"
  , transparency=0.5
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

CASE 2: YOU DO NOT HAVE RAW DATA BUT SYNTHETIC DATA ONLY

... AND YOU WANT TO USE Statistics:-Histogram

Here are the synthetic data which correspond to the raw data R above

>	Classes := [seq(C[k]..C[k+1], k=1..numelems(C)-1)]: SyntheticData := Statistics:-TallyInto(R, Classes): SyntheticData := Matrix(numelems(Classes), 2, (i, j) -> `if`(j=1, lhs(SyntheticData[i]), rhs(SyntheticData[i])));

(1)

If you want to use Statistics:-Histogram you must create hypothetic raw data in a correct way.
The reconstruction below produces raw data that have the same histogram as above.

S := [seq((add(op(SyntheticData[k, 1]))/2)$(SyntheticData[k, 2]), k=1..numelems(Classes))]:

Statistics:-Histogram(
  S
  , frequencyscale=absolute
  , binbounds=C
  , axis[1]=[tickmarks=C]
  , view=[-1..101, default]
  , color="Pink"
  , transparency=0.5
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

CASE 3: YOU DO NOT HAVE RAW DATA BUT SYNTHETIC DATA ONLY

... AND YOU DO NOT NECESSARILLY WANT TO USE Statistics:-Histogram

bars := [
          seq(
            plottools:-rectangle(
              [op(1, SyntheticData[k, 1]), 0]
              , [op(2, SyntheticData[k, 1]), SyntheticData[k, 2]]
              , color="Pink"
              , transparency=0.5
            )
            , k=1..numelems(Classes)
          )
        ]:

plots:-display(
  bars
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

Download Different_options.mw

Be careful when using 2D mode...

Answered: mmcdara 6635

February 18 2024

1 0

It's likely that some typo or hidden character appeared in ode2 (bold red text).
As I never use the document mode with 2D inputs, the attached file displays two types of fonts: so correct your own file the way I did if this annoys you.

restart

>	kernelopts(version)

(1)

(2)

ode1 := diff(f(x), x, x, x)+(1/2*(1+p*(diff(u[e](x), x))/u[e]))*f(x)*(diff(f(x), x, x))+p*(D(u[e]))(1-(diff(f(x), x))^2)/u[e]-p*((diff(f(x), x))*(diff(f(x), x, x))-(diff(f(x), x))*(diff(f(x), x, x))) = 0;

diff(diff(diff(f(x), x), x), x)+(1/2)*(1+p*(diff(u[e](x), x))/u[e])*f(x)*(diff(diff(f(x), x), x))+p*(D(u[e]))(1-(diff(f(x), x))^2)/u[e] = 0

(diff(diff(theta(x), x), x))/Pr+(1/2)*(1+p*(diff(u[e](x), x))/u[e])*f(x)*(diff(theta(x), x))+16*p*a*sigma*L*T[e]^3*(1-theta(x))/(rho*c[p]*u[e]*U[infinity])+U[infinity]^2*p*(diff(f(x), x))*u[e]*(diff(u[e](x), x))/(c[p]*(T[w]-T[e]))-U[infinity]^2*u[e]^2*(diff(diff(f(x), x), x))^2/(c[p]*(T[w]-T[e])) = 0

(3)

bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(8) = 1;

(4)

a1 := [a = 1/2, L = 5, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(T[w]-T[e]) = 10, sigma = 5.67*10^(-8), rho = 1, p = .5];

[a = 1/2, L = 5, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

[a = 1/2, L = 6, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

[a = 1/2, L = 8, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

(5)

>	b1 := eval(ode2, a1); b2 := eval(ode2, a2); b3 := eval(ode2, a3);

(6)

>	c1 := dsolve({bcs1, eval(b1, f(x) = x)}, numeric); c1(0); c2 := dsolve({bcs1, eval(b2, f(x) = x)}, numeric); c2(0); c3 := dsolve({bcs1, eval(b3, f(x) = x)}, numeric); c3(0);

$proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7110816593827365, (3) = 1.4252612944250316, (4) = 2.1450075553641894, (5) = 2.870698607458807, (6) = 3.6016439190865284, (7) = 4.335950433427024, (8) = 5.072381437739985, (9) = 5.809614980858862, (10) = 6.547193597737968, (11) = 7.281757576380688, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .4938484453872593, (2, 1) = .33592640765199755, (2, 2) = .43888481949800257, (3, 1) = .6120099871592314, (3, 2) = .32852292968068963, (4, 1) = .8036666351004551, (4, 2) = .20577920071677938, (5, 1) = .9150410854740033, (5, 2) = .10721924293183477, (6, 1) = .9688640080289379, (6, 2) = 0.4626167885323632e-1, (7, 1) = .990390459943389, (7, 2) = 0.165006856163208e-1, (8, 1) = .9975131515965905, (8, 2) = 0.4861478032432699e-2, (9, 1) = .999462252457013, (9, 2) = 0.11840633668385558e-2, (10, 1) = .9999041167369538, (10, 2) = 0.23851519861008805e-3, (11, 1) = .9999870438515986, (11, 2) = 0.4007921134972896e-4, (12, 1) = 1.0, (12, 2) = 0.58484376487594675e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .4938484453872593, (1, 2) = -0.237191234908843e-1, (2, 1) = .43888481949800257, (2, 2) = -.12498027451545507, (3, 1) = .32852292968068963, (3, 2) = -.1730836385803869, (4, 1) = .20577920071677938, (4, 2) = -.1591461343770899, (5, 1) = .10721924293183477, (5, 2) = -.10974309693499718, (6, 1) = 0.4626167885323632e-1, (6, 2) = -0.59054851436050264e-1, (7, 1) = 0.165006856163208e-1, (7, 2) = -0.25269084092256245e-1, (8, 1) = 0.4861478032432699e-2, (8, 2) = -0.8689730687886674e-2, (9, 1) = 0.11840633668385558e-2, (9, 2) = -0.24203881956517473e-2, (10, 1) = 0.23851519861008805e-3, (10, 2) = -0.5488360802490405e-3, (11, 1) = 0.4007921134972896e-4, (11, 2) = -0.10245379376869548e-3, (12, 1) = 0.58484376487594675e-5, (12, 2) = -0.16375625411613887e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7110816593827365, (3) = 1.4252612944250316, (4) = 2.1450075553641894, (5) = 2.870698607458807, (6) = 3.6016439190865284, (7) = 4.335950433427024, (8) = 5.072381437739985, (9) = 5.809614980858862, (10) = 6.547193597737968, (11) = 7.281757576380688, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.45502255784570956e-7, (2, 1) = 0.9444480543660673e-8, (2, 2) = -0.32701855724739926e-7, (3, 1) = -0.13111776668557004e-6, (3, 2) = 0.7051671015686677e-7, (4, 1) = -0.10085891808028474e-6, (4, 2) = 0.12193667249947745e-6, (5, 1) = 0.13308574074702087e-6, (5, 2) = -0.7660319092840244e-7, (6, 1) = 0.16436590818059128e-6, (6, 2) = -0.18881505625378492e-6, (7, 1) = -0.20116971988762455e-7, (7, 2) = 0.7685426832740135e-9, (8, 1) = -0.9569009446204095e-7, (8, 2) = 0.13118291896195288e-6, (9, 1) = -0.31100640223178834e-7, (9, 2) = 0.4639384647499807e-7, (10, 1) = 0.1584209066120383e-7, (10, 2) = -0.35974084912581454e-7, (11, 1) = 0.11642475203736646e-7, (11, 2) = -0.2764855504318049e-7, (12, 1) = .0, (12, 2) = -0.1925248521887889e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8881505625378492e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8881505625378492e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error end try end proc$

$proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7109194006521007, (3) = 1.4249868973749966, (4) = 2.1446957738695653, (5) = 2.8704181759615643, (6) = 3.6014311331509616, (7) = 4.335815327649273, (8) = 5.07232139321572, (9) = 5.809624985093265, (10) = 6.54727068912088, (11) = 7.281870021183384, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .49766174438535715, (2, 1) = .3375757993634647, (2, 2) = .44007354626176926, (3, 1) = .6138389834751503, (3, 2) = .32815856063407745, (4, 1) = .8049921933517767, (4, 2) = .2049224933509404, (5, 1) = .9157759787447678, (5, 2) = .1064979573437486, (6, 1) = .9691886486674092, (6, 2) = 0.45847467449322866e-1, (7, 1) = .990506554830634, (7, 2) = 0.16320440331941466e-1, (8, 1) = .9975470109277185, (8, 2) = 0.479979457345354e-2, (9, 1) = .9994703351484253, (9, 2) = 0.1167147574584417e-2, (10, 1) = .9999056808640605, (10, 2) = 0.23475927713373983e-3, (11, 1) = .9999872681820644, (11, 2) = 0.3939868164434801e-4, (12, 1) = 1.0, (12, 2) = 0.5746972025352357e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .49766174438535715, (1, 2) = -0.279454400456164e-1, (2, 1) = .44007354626176926, (2, 2) = -.12801162336374522, (3, 1) = .32815856063407745, (3, 2) = -.1744590166938988, (4, 1) = .2049224933509404, (4, 2) = -.15927332083416956, (5, 1) = .1064979573437486, (5, 2) = -.10934646266608739, (6, 1) = 0.45847467449322866e-1, (6, 2) = -0.58651810580892944e-1, (7, 1) = 0.16320440331941466e-1, (7, 2) = -0.25032143866202835e-1, (8, 1) = 0.479979457345354e-2, (8, 2) = -0.8589685100784845e-2, (9, 1) = 0.1167147574584417e-2, (9, 2) = -0.23880431153508012e-2, (10, 1) = 0.23475927713373983e-3, (10, 2) = -0.5405971765590906e-3, (11, 1) = 0.3939868164434801e-4, (11, 2) = -0.10076942378371007e-3, (12, 1) = 0.5746972025352357e-5, (12, 2) = -0.16091521666157516e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7109194006521007, (3) = 1.4249868973749966, (4) = 2.1446957738695653, (5) = 2.8704181759615643, (6) = 3.6014311331509616, (7) = 4.335815327649273, (8) = 5.07232139321572, (9) = 5.809624985093265, (10) = 6.54727068912088, (11) = 7.281870021183384, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.4573949665850558e-7, (2, 1) = 0.6242874960547203e-8, (2, 2) = -0.3057749989248847e-7, (3, 1) = -0.12966272475504071e-6, (3, 2) = 0.7217465704512478e-7, (4, 1) = -0.940390491554525e-7, (4, 2) = 0.1178112154625901e-6, (5, 1) = 0.1357038240627767e-6, (5, 2) = -0.8040712939043494e-7, (6, 1) = 0.1607319279832156e-6, (6, 2) = -0.18568759788316905e-6, (7, 1) = -0.2282862282805592e-7, (7, 2) = 0.4786547870156748e-8, (8, 1) = -0.94843869899127e-7, (8, 2) = 0.13055152558459717e-6, (9, 1) = -0.2970037963750006e-7, (9, 2) = 0.4434338125388376e-7, (10, 1) = 0.16200082159941674e-7, (10, 2) = -0.36390277880281064e-7, (11, 1) = 0.11561105767115588e-7, (11, 2) = -0.27217683490459804e-7, (12, 1) = .0, (12, 2) = -0.16987712182963502e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8568759788316905e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8568759788316905e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error end try end proc$

$proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7105918488760237, (3) = 1.4244328611291581, (4) = 2.1440659423124293, (5) = 2.8698512980998347, (6) = 3.6010004098896413, (7) = 4.335541208202285, (8) = 5.072198845420454, (9) = 5.809644194733494, (10) = 6.5474258600340365, (11) = 7.282096976039094, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .5052552537006175, (2, 1) = .3408462755204489, (2, 2) = .442422215142914, (3, 1) = .6174560224034431, (3, 2) = .32742681953867314, (4, 1) = .8076074937875736, (4, 2) = .2032242438976165, (5, 1) = .9172227818187347, (5, 2) = .10507373726470326, (6, 1) = .9698264569016594, (6, 2) = 0.45031941250935854e-1, (7, 1) = .9907341917436213, (7, 2) = 0.15966437963481905e-1, (8, 1) = .9976132778134335, (8, 2) = 0.4678915550247515e-2, (9, 1) = .9994861265216441, (9, 2) = 0.11340645752511988e-2, (10, 1) = .9999087319608427, (10, 2) = 0.22742703929550063e-3, (11, 1) = .9999877052438856, (11, 2) = 0.38072228239622816e-4, (12, 1) = 1.0, (12, 2) = 0.5549311882673039e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .5052552537006175, (1, 2) = -0.364325736980702e-1, (2, 1) = .442422215142914, (2, 2) = -.13404823355722884, (3, 1) = .32742681953867314, (3, 2) = -.17717619408563023, (4, 1) = .2032242438976165, (4, 2) = -.15951351711336537, (5, 1) = .10507373726470326, (5, 2) = -.10855688752010623, (6, 1) = 0.45031941250935854e-1, (6, 2) = -0.5785531343159658e-1, (7, 1) = 0.15966437963481905e-1, (7, 2) = -0.2456567964351268e-1, (8, 1) = 0.4678915550247515e-2, (8, 2) = -0.8393290947591843e-2, (9, 1) = 0.11340645752511988e-2, (9, 2) = -0.2324700819301251e-2, (10, 1) = 0.22742703929550063e-3, (10, 2) = -0.5244967168305387e-3, (11, 1) = 0.38072228239622816e-4, (11, 2) = -0.9748390992631722e-4, (12, 1) = 0.5549311882673039e-5, (12, 2) = -0.15538073266825603e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7105918488760237, (3) = 1.4244328611291581, (4) = 2.1440659423124293, (5) = 2.8698512980998347, (6) = 3.6010004098896413, (7) = 4.335541208202285, (8) = 5.072198845420454, (9) = 5.809644194733494, (10) = 6.5474258600340365, (11) = 7.282096976039094, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.4621216507042107e-7, (2, 1) = 0.18871763805575788e-9, (2, 2) = -0.26440976380442328e-7, (3, 1) = -0.12642870795085007e-6, (3, 2) = 0.7507398457807075e-7, (4, 1) = -0.8067853301949186e-7, (4, 2) = 0.10949148171767816e-6, (5, 1) = 0.1404802397309745e-6, (5, 2) = -0.8760075826073426e-7, (6, 1) = 0.15342379606531546e-6, (6, 2) = -0.1792563287988756e-6, (7, 1) = -0.28042403695503607e-7, (7, 2) = 0.12597322677227419e-7, (8, 1) = -0.9307097617997284e-7, (8, 2) = 0.12914164433307568e-6, (9, 1) = -0.26946194672378708e-7, (9, 2) = 0.4029811073645486e-7, (10, 1) = 0.1687782559135978e-7, (10, 2) = -0.37161040827639946e-7, (11, 1) = 0.11392020809816244e-7, (11, 2) = -0.2635259175010777e-7, (12, 1) = .0, (12, 2) = -0.12553181464882845e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.792563287988756e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.792563287988756e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error end try end proc$

(7)

p1 := odeplot(c1, [x, theta(x)], 0 .. 8, colour = green):

Download PGVTN_mmcdara.mw

PS 1: using PDETools:-declare doesn't ease the task of finding where the error is and you must be doubly careful using it.
PS 2: be consistent and declare u[e] too if you consider it as a function of x (which u[e]' seems to mean, even if u[e] is later defined as a constant function).

Use orthogonality condition of Hermite p...

Answered: mmcdara 6635

February 17 2024

2 0

As explained in the attached file your code does not work with my Maple 2015?
So I copy-pasted your expression of eq5 and used it as the starting point to get eq6.

Your code is in grey, mine on bold brown; the successive transformations of F (the integral term in your eq6) enableobtaoining a standard form from which it clearly appears that the orthogonality condition of Hermite polynomials (Maple uses their physical form) can be used to simplify the expression of F:

Final remark: as _C2 intervenes through its square, there are two real and opposite values of _C2 such that F(n)=1.

restart;

>	kernelopts(version)

(1)

Here is your file NOT EXECUTED in Maple 2015 for PDEtools:-dchange(...) generates an error with Maple 2015

>	with(PDEtools): #with(Maplets[Elements]):

>	eq0:= E1 = (n+1/2)_h; param := m1 =_hsigma^2;

(2)

>	os := expand(isolate((-_h^2/(2m1))diff(psi(x),x,x)+(m1x^2/2)psi(x)-E1*psi(x)=0,diff(psi(x),x,x)));

diff(diff(psi(x), x), x) = m1^2*x^2*psi(x)/_h^2-2*m1*E1*psi(x)/_h^2

(3)

eq1 := PDEtools:-dchange([x=sqrt(_h/m1)*zeta,psi(x)=f(zeta)*exp(-zeta^2/2)],os,[zeta,f(zeta)],params={m1,_h,E1});
eq2 := expand(isolate(eq1,diff(f(zeta),zeta$2)));
eq3 := simplify(eval(remove(has,convert(rhs(dsolve(eq2,f(zeta))),Hermite),KummerM),eq0)*exp(-zeta^2/2),radical) assuming zeta::positive;
eq4 := subs(param,simplify(subs(zeta=sqrt(_h/m1)*x,eq3))) ;
eq5 := unapply(eq4,n,x,sigma);
eq6 := solve(int(eq5(n,x,sigma)^2, x=-infinity...infinity)=1,_C2,useassumptions) assuming _C2 :: real, sigma::positive,n::posint;

((diff(diff(f(zeta), zeta), zeta))*exp(-(1/2)*zeta^2)-2*(diff(f(zeta), zeta))*zeta*exp(-(1/2)*zeta^2)-f(zeta)*exp(-(1/2)*zeta^2)+f(zeta)*zeta^2*exp(-(1/2)*zeta^2))*m1/_h = m1*zeta^2*f(zeta)*exp(-(1/2)*zeta^2)/_h-2*m1*E1*f(zeta)*exp(-(1/2)*zeta^2)/_h^2

diff(diff(f(zeta), zeta), zeta) = -2*E1*f(zeta)/_h+2*(diff(f(zeta), zeta))*zeta+f(zeta)

_C2*HermiteH(n, zeta)*exp(-(1/2)*zeta^2)/2^n

_C2*HermiteH(n, (1/sigma^2)^(1/2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2)

proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, (1/sigma^2)^(1/2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

(4)

From now on: I copy-paste the expression of "your" (that I cannot obtain with Maple 2015)

>	eq5 := proc (n, x, sigma) options operator, arrow; _C2HermiteH(n, sqrt(1/sigma^2)x)2^(-n)exp(-(1/2)*x^2/sigma^2) end proc

proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, sqrt(1/sigma^2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

(5)

Here is the integral contained in your relation "eq6"

After some elementary algebra I get the last relation above

>	F := Int(eq5(n,x,sigma)^2, x=-infinity...infinity) assuming n::posint; F := combine(F, exp); F := eval(F, sqrt(1/sigma^2) = 1/sigma); F := simplify(IntegrationTools:-Expand(F)); F := IntegrationTools:-Change(F, x/sigma=z, z) assuming sigma > 0; F := simplify(IntegrationTools:-Expand(F));

Int(_C2^2*HermiteH(n, (1/sigma^2)^(1/2)*x)^2*(2^(-n))^2*(exp(-(1/2)*x^2/sigma^2))^2, x = -infinity .. infinity)

Int(_C2^2*HermiteH(n, (1/sigma^2)^(1/2)*x)^2*(2^(-n))^2*exp(-x^2/sigma^2), x = -infinity .. infinity)

Int(_C2^2*HermiteH(n, x/sigma)^2*(2^(-n))^2*exp(-x^2/sigma^2), x = -infinity .. infinity)

_C2^2*(Int(HermiteH(n, x/sigma)^2*exp(-x^2/sigma^2), x = -infinity .. infinity))*4^(-n)

_C2^2*(Int(HermiteH(n, z)^2*exp(-z^2)*sigma, z = -infinity .. infinity))*4^(-n)

_C2^2*sigma*(Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity))*4^(-n)

(6)

From orthogonality of Hermite polynomials (Physical form used in Maple)
(see Wiki)

>	select(has, indets(F, function), Int)[]; OrthogonalityCondition := % = (2^nn!sqrt(Pi))

Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity)

Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity) = 2^n*factorial(n)*Pi^(1/2)

(7)

>	F := unapply(eval(F, OrthogonalityCondition), n)

(8)

Thus the expression of _C2 for any n

>	eq6 := unapply( [solve(F(n)=1, _C2)], n)

proc (n) options operator, arrow; [1/(factorial(n)*Pi^(1/2)*2^(-n)*sigma)^(1/2), -1/(factorial(n)*Pi^(1/2)*2^(-n)*sigma)^(1/2)] end proc

(9)

Download SolnforOS_dsolve_mmcdara.mw

The two derivatives give exactly the sam...

Answered: mmcdara 6635

February 14 2024

1 2

The two derivatives give exactly the same thing a.

Nevertheless here is a way to do what you want derivatives.mw

Personally...

Answered: mmcdara 6635

February 10 2024

1 3

I_would_do_that.mw

I declared prnt and clr as optional parameters (first one is boolean initialized to true, second one is string initialized to blue), I introduced a control of the color you pass, and I simplified your code a little bit.
Here are the modifications (complete code in the attached file)

spread2:=proc(
               ....
               {clr::string:="Blue"}, 
               {prnt::boolean:=true}
             )
           ....
           local col := substring(StringTools:-Capitalize(clr), 1..1);
           local COL := `if`(col="R", "Red", `if`(col="G", "Green", "Blue"))  
           ....
           if `not`(member(col, {"R", "G", "B"})) then
             error cat("allowed colors are red, green or blue, received ", clr)
           end if:

           if prnt then
             print(cat("Spread 2 [x,y] Points/Vectors wrt origin ",  COL));
           end if:
                 
           if col = "B" then
               return ....

           elif col="G" then
               return ....

           elif col="R" then
               return ....
          end if;
          end proc:

More controls could be necessary for a safer code (what if the denominators in the green and red cases are null?)

A way...

Answered: mmcdara 6635

February 08 2024

2 0

Download arctan.mw

Unfortunately I gave only Maple 2015 and...

Answered: mmcdara 6635

February 03 2024

1 1

Unfortunately I gave only Maple 2015 and WeibullPlot is not avaliable.

Look to this file and tell me if the plot I draw corresponds or not to what you would like to get?
I'm afraid I did not corerctly understood what is the issue you refer to.
Here axis 1 is given 4 subticks and axis 2 3 subricks, you can see that the plot is correct.

WeibullPlot(Vector[row](%id = 18446744078086407582), reference = false, style = line, gridlines = true, size = [800, 500], axis = [tickmarks = [default, subticks = 4], color = darkgreen])

(1)

>	kernelopts(version)

(2)

F := unapply(CDF(Weibull(1, 0.6), z), z):

AA := sort(AA):
N := numelems(AA):

interface(warnlevel=0):
WB := map(n -> [AA[n], solve(F(z)=n/N)], [$1..N-1]):

plots:-display(
  plot(
    WB
    , color="NavyBlue"
    , gridlines=true
    , axis[1]=[mode=log, color="DarkGreen", tickmarks=[subticks=3]]
    , axis[2]=[mode=log, color="DarkGreen", tickmarks=[subticks=4]]
    , size=[700, 500]
    , axes=boxed
  )
  , plot(Mean(Weibull(1, 0.6)), z=(min..max)(AA), color=red, linestyle=3)
  , plot([[Mean(AA), WB[1][2]], [Mean(AA), WB[-1][2]]], z=(min..max)(AA), color=red, linestyle=3)
)

>	Sample_data := log[10]~(map2(op, 1, WB)): Reference_data := log[10]~(map2(op, 2, WB)): FIT := LinearFit([1, x], Reference_data, Sample_data, [x], output=solutionmodule): print~(FIT:-Results()):

"parametervalues" = (Vector(2, {(1) = 0.566005759094454e-1, (2) = .996156636676824}))

"parametervector" = (Vector(2, {(1) = 0.566005759094454e-1, (2) = .996156636676824}))

"confidenceintervals" = (Vector(2, {(1) = 0.413122294090191e-1 .. 0.718889224098718e-1, (2) = .980270809958746 .. 1.01204246339490}))

$"residuals" = (Vector(4, {(1) = ` 1 .. 99 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}))$

$"leverages" = (Vector(4, {(1) = ` 1 .. 99 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}))$

"variancecovariancematrix" = (Matrix(2, 2, {(1, 1) = 0.593363937315893e-4, (1, 2) = 0.259631230458869e-4, (2, 1) = 0.259631230458869e-4, (2, 2) = 0.640648414147044e-4}))

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

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

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

"AtkinsonTstatistic" = Vector[column](%id = 18446744078229089086)

(3)

p := plots:-display(
  plot([seq([Reference_data[n], Sample_data[n]], n=1..numelems(Sample_data))], color="NavyBlue")
  , plot(
      FIT:-Results("leastsquaresfunction"), x=(min..max)(Reference_data), color="NavyBlue", linestyle=3
    )
  , gridlines=true
):

p := plottools:-transform((x, y) -> [y, x])(p):

plots:-display(
  p
  , axis[1]=[tickmarks=[seq(k=10.0^k, k=-3..1)], color="DarkGreen"]
  , axis[2]=[tickmarks=[seq(k=10.0^k, k=-3..1)], color="DarkGreen"]
  , view=[(min..max)(Sample_data), (min..max)(Reference_data)]
  , size=[700, 500]
  , axes=boxed
)

Download WeibullPlot_2015.mw

The limit is not a number...

Answered: mmcdara 6635

February 02 2024

0 2

but an interval.

Just do this to see what happens for a given y as x goes to 0.

limit(f(x, y), x=0)
                            /1\             /1\
                   -1 + sin| - | .. 1 + sin| - |
                            \y/             \y/

Maybe doing this will help you understand what happens

limit(f(1/x, 1/y), x=+infinity)
                   -1 + sin(y) .. 1 + sin(y)
limit(sin(x), x=+infinity)
                            -1 .. 1

You can also compute a directional limit

limit(f(x, lambda*x), x=0)
                            -2 .. 2

In the attached file you will find an exploration of f(x,y) as you get closer and closer to the point (0, 0) (the value h of the slider means the neighborhood of (0, 0) has size 10^-h).
limit.mw

A quick slution...

Answered: mmcdara 6635

February 02 2024

0 2

A quick workaround is

plots:-implicitplot(g=1e-8, x = 0 .. 80, y = 0 .. 60, scaling = constrained, grid=[400, 400])

Nevertheless do not gorget that g=0 on the domain 33 <= x <= 58, 0, <= y <= 49.99226668
gnull.mw

One way...

Answered: mmcdara 6635

January 22 2024

1 1

(works only is the pattern to search is numeric, if not adaptations have to be done)

restart

>	#M3 := Matrix(4, 3, [34, 67, 1, 35, 80, 1, 45, 78, 2, 56, 99, 2]): M3 := LinearAlgebra:-RandomMatrix(10, 5, generator=0..10)

(1)

LookUp := proc(M::Matrix, cr, P, x)
  local S:
  if cr='column' then
    S := select((n -> is(M[n, P]=x)), [$1..numelems(M[..,1])]);
    if S <> [] then
      return M[S, ..];
    else
      error cat("Column ", P, " does not contain ", convert(x, string))
    end if:
  elif cr='row' then
    S := select((n -> is(M[P, n]=x)), [$1..numelems(M[1, ..])]);
    if S <> [] then
      return M[.., S];
    else
      error cat("Row ", P, " does not contain ", convert(x, string))
    end if:
  end if:
end proc:

>	# Submatrix made of rows whose column "P" contains value "x" LookUp(M3, 'column', 5, 4); LookUp(M3, 'column', 2, 1);

Matrix([[0, 2, 6, 5, 4], [3, 10, 4, 9, 4], [7, 6, 6, 2, 4]])

(2)

>	# Submatrix made of columns whose row "P" contains value "x" LookUp(M3, 'row', 2, 10); LookUp(M3, 'row', 1, 4);

Matrix([[2, 6], [10, 10], [5, 5], [6, 1], [5, 3], [10, 4], [1, 6], [2, 5], [0, 3], [6, 6]])

(3)

>	type(LookUp(M3, 'row', 2, 31), numeric)

Error, (in LookUp) Row 2 does not contain 31

Download LookUp.mw

Note that I did not pay attention to performance and that the LookUp function may be inefficient for large matrices.

In case you would like to remove the columns/rows which contain the pattern to search, here is variant of the previous file LookUp_variant.mw

If I understand correctly your problem.....

Answered: mmcdara 6635

January 18 2024

2 6

... here is a solution I wonder what you mean by "The result will be a list of 56 models with 5 monomials."?)

It is efficient while the lengths of the two lists are not much larger:

restart:
with(combinat):
L1 := [x^2*y*alpha[1, 11], x*z^2*alpha[2, 15], y^2*z*alpha[3, 17] + y*z*alpha[3, 8]]:
         
L2 := [alpha[i, 0], alpha[i, 1]*x, alpha[i, 2]*y, alpha[i, 3]*z, alpha[i, 4]*x^2, alpha[i, 5]*y*x, alpha[i, 6]*z*x, alpha[i, 7]*y^2, alpha[i, 8]*z*y, alpha[i, 9]*z^2, alpha[i, 10]*x^3, alpha[i, 11]*y*x^2, alpha[i, 12]*z*x^2, alpha[i, 13]*y^2*x, alpha[i, 14]*z*y*x, alpha[i, 15]*z^2*x, alpha[i, 16]*y^3, alpha[i, 17]*z*y^2, alpha[i, 18]*z^2*y, alpha[i, 19]*z^3]:

#-----------------------------------------------------------------

# Build the cartesian product of the two lists

T:=cartprod([L1, L2]):
L1L2 := NULL:
while not T[finished] do 
  L1L2 := L1L2, T[nextvalue]() 
end do:
L1L2 := [L1L2]:

# Keep only the elements which do not have repeated alphas with the same second index.

map(u -> map2(op, 2, [select(type, indets(u), `indexed`)[]]), L1L2):
Keep := zip((u, v) -> if `not`(member(v[-1], {v[1..-2][]})) then add(u) end if, L1L2, %):

print~ (Keep): # to get Keep displayed one element per line

numelems(Keep)
                               56

A_solution.mw

If your matrix has the structure you giv...

Answered: mmcdara 6635

January 04 2024

3 1

If your matrix has the structure you give, that is +/- a given quantity at some positions and 0 elsewhere, here is a way=

>	kernelopts(version)

(1)

>	b := Matrix(3, 6, [[-1/2, 0, 1/2, 0, 0, 0], [0, 0, 0, -1/2, 0, 1/2], [0, -1/2, -1/2, 1/2, 1/2, 0]]);

b := Matrix(3, 6, {(1, 1) = -1/2, (1, 2) = 0, (1, 3) = 1/2, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = -1/2, (2, 5) = 0, (2, 6) = 1/2, (3, 1) = 0, (3, 2) = -1/2, (3, 3) = -1/2, (3, 4) = 1/2, (3, 5) = 1/2, (3, 6) = 0})

(2)

>	opb := op(b)[3];

(3)

>	ropb := rhs~(opb)

(4)

>	b1 := ropb[1] * ``(b / ropb[1])

b1 := -Typesetting[delayDotProduct](1/2, Matrix(3, 6, {(1, 1) = 1, (1, 2) = 0, (1, 3) = -1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (2, 5) = 0, (2, 6) = -1, (3, 1) = 0, (3, 2) = 1, (3, 3) = 1, (3, 4) = -1, (3, 5) = -1, (3, 6) = 0}), true)

(5)

>	b2 := ropb[2] * ``(b / ropb[2])

Typesetting[delayDotProduct](1/2, Matrix(3, 6, {(1, 1) = -1, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = -1, (2, 5) = 0, (2, 6) = 1, (3, 1) = 0, (3, 2) = -1, (3, 3) = -1, (3, 4) = 1, (3, 5) = 1, (3, 6) = 0}), true)

(6)

Download For_instance.mw

For more complex matrices there are some ambiguities about what is the "common factor" to put outside of the matrix: For_instance_2.mw

This way...

Answered: mmcdara 6635

January 03 2024

3 0

restart

# list of small greek letters
seq(cat(`#mo("&#`, i, `;")`), i=945..969)



# A procedure to generate a random word of L letters
RandomGreekWord := proc(L)
  local r := rand(945..969);
  cat(
       `#mrow(`, 
       seq(cat(`mo("&#`, r(), `;"),`), i=1..L-1),
       cat(`mo("&#`, r(), `;"))`)
  )
end proc:

# an example
RandomGreekWord(10)

Download Greek.mw

Capital greek letters range from 913 to 937.
So all greek letters can be displayed by replacing

i=945..969

i in GreekIndices

where

GreekIndices := [$913..937, $945..969]

Several possibilities...

Answered: mmcdara 6635

December 28 2023

2 1

I didn(t clearly understood if you had only 3 set of parameters (your second question) or 27 (3x3x3).
Nevertheless the two cases are adressed in the attached file.

As you load the Statistics package I guessed you intended to use it somewhere. So you will find in this same file a very short illustration of what you could do.

Solving an ODE in maple is straigthforward

(1)

(2)

(3)

(4)

ODE1 := diff(P(x), x) = a1*(P__b1-P(x))/(1+lambda1(P(x)))^(2/3)

(5)

(6)

(7)

(8)

First question how can I get only the value of P(x) at LLast?

>	eval(P(x), sol1(L__Last)); select(has, sol1(L__Last), P(x)); select(has, sol1(L__Last), P(x))[];

(9)

Second question how can I solve this ODE for multiple input and finally obtain a vector with all results?

>	# my prefered way is this one restart: ODE1 := diff(P(x), x) = a1(P__b1-P(x))/(1+lambda1(P(x)))^(2/3); params := convert(indets(ODE1, name) minus {x}, list); sol1 := dsolve({ODE1, P(0) = 40.42}, numeric, parameters=params);

diff(P(x), x) = a1*(P__b1-P(x))/(1+lambda1*P(x))^(2/3)

(10)

>	data_Pb1 := [284, 283, 352]: data_a1 := [0.150, 0.152, 0.145]: data_lam := [38.57, 50.22, 39.83]: data := [ data_Pb1, data_a1, data_lam ]

(11)

>	# To get the 3x3x3=27 different combinations of the three parameters: use combinat in C := cartprod(data): combs := NULL: while not C[finished] do combs := combs, C[nextvalue]() end do: end use: combs := [combs]:

>	# To get the 27 corrresponding values: L__Last := 10: RES := table([]): # a suggestion for c in combs do sol1(parameters = c): RES[c] := eval(P(x), sol1(L__Last)); end do: eval(RES):

# Visualizations.
#
# For instance, to plot P(L__last) as a "function" of P__b1 and a1
# "parameterized" by λ1

# Indices of table RES such that λ1 = 38.57
Indices_1 := select((i -> is(op(3, i)=data_lam[1])), [indices(RES, nolist)]):

# Indices of table RES such that λ1 = 50.22
Indices_2 := select((i -> is(op(3, i)=data_lam[2])), [indices(RES, nolist)]):

# Indices of table RES such that λ1 = 39.83
Indices_3 := select((i -> is(op(3, i)=data_lam[3])), [indices(RES, nolist)]):

# Results for λ1 = 39.83
N_1 := numelems(Indices_1):
pts_1 := Matrix(N_1, 3, (i, j) -> `if`(j=3, RES[Indices_1[i]], Indices_1[i][j])):

# Results for λ1 = 39.83
pts_2 := Matrix(N_1, 3, (i, j) -> `if`(j=3, RES[Indices_2[i]], Indices_2[i][j])):
# Results for λ1 = 39.83
pts_3 := Matrix(N_1, 3, (i, j) -> `if`(j=3, RES[Indices_3[i]], Indices_3[i][j])):

>	plots:-display( Statistics:-ScatterPlot3D(pts_1, symbol=solidcircle, symbolsize=25, color=red), Statistics:-ScatterPlot3D(pts_2, symbol=solidcircle, symbolsize=25, color=green), Statistics:-ScatterPlot3D(pts_3, symbol=solidcircle, symbolsize=25, color=blue) );

ResIndices := [indices(RES, nolist)]:
ResMatrix := convert(ResIndices, Matrix):
ResMatrix := `<|>`(ResMatrix, Vector(numelems(combs), i -> RES[ResIndices[i]]));

fit := unapply(Statistics:-LinearFit([1, u, v, w], ResMatrix, [u, v, w]), [u, v, w]):

fit(params[]);

plots:-display(
  Statistics:-ScatterPlot3D(pts_1, symbol=solidcircle, symbolsize=25, color=red),
  Statistics:-ScatterPlot3D(pts_2, symbol=solidcircle, symbolsize=25, color=green),
  Statistics:-ScatterPlot3D(pts_3, symbol=solidcircle, symbolsize=25, color=blue),

  plot3d(fit(P__b1, a1, data_lam[1]), P__b1=(min..max)(data_Pb1), a1=(min..max)(data_a1), style=wireframe, color=red),
  plot3d(fit(P__b1, a1, data_lam[2]), P__b1=(min..max)(data_Pb1), a1=(min..max)(data_a1), style=wireframe, color=green),

  plot3d(fit(P__b1, a1, data_lam[3]), P__b1=(min..max)(data_Pb1), a1=(min..max)(data_a1), style=wireframe, color=blue)
);

$ResMatrix := Vector(4, {(1) = ` 27 x 4 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})$

# In case you have just 3 configurations:

ThreeRes := convert([data_Pb1, data_a1, data_lam], Matrix)^+:
ThreeConf := convert(ThreeRes, listlist):

ResVector := Vector(numelems(ThreeConf)):
for i from 1 to numelems(ThreeConf) do
sol1(parameters = ThreeConf[i]):
ResVector[i] := eval(P(x), sol1(L__Last));
end do:

VectorTitle := Vector[row]([params[], P(L__Last)]):

ThreeRes := `<,>`(VectorTitle, < ThreeRes | ResVector >)

ThreeRes := Matrix(4, 4, {(1, 1) = `#msub(mi("P"),mi("b1"))`, (1, 2) = a1, (1, 3) = lambda1, (1, 4) = P(10), (2, 1) = 284, (2, 2) = .150, (2, 3) = 38.57, (2, 4) = 43.0642521778387, (3, 1) = 283, (3, 2) = .152, (3, 3) = 50.22, (3, 4) = 42.6670938150578, (4, 1) = 352, (4, 2) = .145, (4, 3) = 39.83, (4, 4) = 43.6076477820810})

(12)

Download SeveralPossibilities.mw

You need to cheat...

Answered: mmcdara 6635

December 27 2023

4 2

@sursumCorda

Basically you have four elements and each is associated to a frame: to satisfy your requirement the red line must be set behind the blue rectangle, the blue rectangle behind the blue line, finally the blue line behind the red rectangle... which implies the red line is behind the red rectangle.
So it is impossible to do what you asked without cheating

Here are two ways to cheat (is the first way what you name "breaking segments"?)

restart

>	with(ColorTools): with(plottools):

>	lr := rectangle([2, 0], [3, 6], 'color' = Color([.6, .7, .9])):

>	rr := rectangle([4, 0], [5, 6], 'color' = Color([.95, .7, .6])):

>	h := 0.1: al := [[1, 0], [6, 5]]: al := [al[], (al[[2, 1]] +~ [[0, h]$2])[]]: pl := PLOT(POLYGONS(al, COLOR(RGB, HexToRGB24("#0072BD")[]), STYLE(PATCHNOGRID))): ar := [[1, 1], [6, 6]]: ar := [ar[], (ar[[2, 1]] +~ [[0, h]$2])[]]: pr := PLOT(POLYGONS(ar, COLOR(RGB, 1, 0, 0), STYLE(PATCHNOGRID))):

>	h := 0.1: al := [[1, 0], [6, 5]]: al := [al[], (al[[2, 1]] +~ [[0, h]$2])[]]: lp := PLOT(POLYGONS(al, COLOR(RGB, HexToRGB24("#0072BD")[]), STYLE(PATCHNOGRID))): ar := [[1, 1], [6, 6]]: ar := [ar[], (ar[[2, 1]] +~ [[0, 2*h]$2])[]]: rp := PLOT(POLYGONS(ar, COLOR(RGB, 1, 0, 0), STYLE(PATCHNOGRID))):

>	half_rp := transform((x, y) -> [(x+1)/2, y/2])(rp): plots:-display( translate(half_rp, 5/2, 5/2), rr, lp, lr, rr, half_rp )

While_cheating.mw

Another solution is to build 3D elements in different planes and project the result onto the appropriate plane. For instance:

bl := plot([[1, 0], [6, 5]], color="#0072BD", thickness=3): 

rl := PLOT3D(
        CURVES(
          [[1, 1, -1], [7/2, 7/2, -1], [7/2, 7/2, 1],  [6, 6, 1]]
          , COLOR(RGB, 1, 0, 0)
          , THICKNESS(3)
        )
      ): 

lr3D := transform((x, y) -> [x, y, -1/2])(lr):
rr3D := transform((x, y) -> [x, y, 1/2])(rr):
bl3D := transform((x, y) -> [x, y, 0])(bl):

plots:-display(lr3D, rr3D, bl3D, rl, lightmodel=none, orientation=[-90, 0, 0], labels=[x, y, z])

5 6 7 8 9 10 11

Last

Page 7 of 57

Share via:

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

mmcdara

6635 Reputation

18 Badges

MaplePrimes Activity

These are answers submitted by mmcdara

I hope this will enlighten you enough to...

Be careful when using 2D mode...

Use orthogonality condition of Hermite p...

The two derivatives give exactly the sam...

Personally...

A way...

Unfortunately I gave only Maple 2015 and...

The limit is not a number...

A quick slution...

One way...

If I understand correctly your problem.....

If your matrix has the structure you giv...

This way...

Several possibilities...

You need to cheat...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again