## 6279 Reputation

7 years, 325 days

## One way...

(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);
 (2)
 > # Submatrix made of columns whose row "P" contains value "x" LookUp(M3, 'row', 2, 10); LookUp(M3, 'row', 1, 4);
 (3)
 > type(LookUp(M3, 'row', 2, 31), numeric)
 >

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.....

... 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...

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]]);
 (2)
 > opb := op(b)[3];
 (3)
 > ropb := rhs~(opb)
 (4)
 > b1 := ropb[1] * ``(b / ropb[1])
 (5)
 > b2 := ropb[2] * ``(b / ropb[2])
 >
 (6)
 >

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...

```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)
```

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

`i=945..969`

by

`i in GreekIndices`

where

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

## Several possibilities...

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)

 (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);
 (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) );
 > # 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 >)
 (12)

## You need to cheat...

@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])

```

## 0.*I is not the only problem...

Even if you get rid of 0.*I you will face some other problems (for instance eq1 and eq2 are complex, so my STEP 2 in the attached file)

Look to the attached file to see how I fixed them successively (yellow highlighted text)  while leaving to you the last one (uou are the one who can fix it).

(I did the things quite rapidly and didn't spend time to arrange your code)

## Use Fourier transform's properties...

For the first one you have to use the translation and scaling properties of FT/FFT/DFT transforms.
For the second one you cannot keep a formal N and you must use the scaling property above.

For whatever reason the content of the attached file can't be displayed

FT_mmcdara.mw

A few details.mw here

## One way to get matrix forms...

You could do this

 >
 >
 >
 >
 >
 >
 >
 >
 > kernelopts(version);  # Newer versions ebable a better rendering while avoiding &.
 (1)
 > # A1 A1_rel   := map(a -> isolate(a, select(has, indets(a, function), diff)[]), [entries(A1, nolist)]): Unknowns := lhs~(%): mat, vec := LinearAlgebra:-GenerateMatrix(A1_rel, Unknowns): # This mat &. vec = < Unknowns >: # or this mat . LinearAlgebra:-DiagonalMatrix(vec) = < Unknowns >;
 (2)
 > # A3 Unknowns := convert(select(has, indets(A3), g), list); mat, vec := LinearAlgebra:-GenerateMatrix([entries(A3, nolist)], Unknowns): 'A3' = mat &. < Unknowns >;
 (3)
 > # Verify that lhs = rhs eval(%, `&.` = `.`): simplify([entries(lhs(%), nolist)] -~ [entries(rhs(%), nolist)]);
 (4)
 >

## Would that suit you?...

```restart
with(Units):
UnitFactors_1 := proc(F)
local U, M, u, d:
U := [op(parse(convert(op(F), string)))];
M := 1:
for u in U do
if u::`^` then
d := ldegree(u):
if d < 0 then
M := M / Unit(op(1, u))^(-d):
else
M := M * Unit(op(1, u))^(d):
end if:
else
M := M*Unit(u):
end if;
end do:
M;
end proc:

Newton := Unit('kg'*'m'/'s'^2):

lprint(Newton = UnitFactors_1(Newton));
Units:-Unit(kg*m/s^2) = Units:-Unit(kg)*Units:-Unit(m)/Units:-Unit(s)^2

```

Based upon

`GetUnit(some_compound_unit, conversion)`

this file gives some hints to answer your second question but it is Far_from_perfect.mw

## What if you do this?...

```r1 := -1 <= x and x <= 0:
r2 := 0 <= x and x <= 1:
x  := solve({convert(r1, 'relation') or convert(r2, 'relation')})
{-1 <= x, x <= 1}

`and`(s[])
-1 <= x and x <= 1

```

## A matter of display?...

In Excel the decimal separator is only a matter of display.
If you only want to display Maple integers with a comma (and not do calculus, see @Thomas Richard's answer) you can do this

 > restart
 > # Comma Separated Integer CSI := proc(x)   StringTools:-Substitute(convert(x, string), ".", ","):   nprintf("%s", %) end proc:
 > rand(-10.0 .. 10.0)(); CSI(%);
 (1)
 > # two equivalent forms V  := LinearAlgebra:-RandomVector(5, generator=-100.0 .. 100.0): VC := Vector(5, i -> CSI(V[i])): V, VC , `or`, CSI~(V)
 (2)
 > # With a control of the number of digits CSI := proc(x, n)   local s, sl, sr:   uses StringTools:   s := convert(x, string):   if x::integer then     s := nprintf("%s", s)   else     s := Split(s, "."):     nprintf("%s,%s", %[1], substring(%[2], 1..n)):   end if: end proc:
 > rand(-10.0 .. 10.0)(); CSI( %, 8); CSI(%%, 3);
 (3)
 >

## Your code is full of errors...

Read carefully the attached file to understand where your errors are and how to correct them.
red texts explain where the previous errors come from.

 > restart

Successive errors and how to get rid of them

 > Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j));       p := plot(finite_element_sigma[i](x), x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := display(Graph, p)     end if   end do: end do:
 > # In the outer loop N is undefined N := 10:  # it's up to you to provide the good value Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j));       p := plot(finite_element_sigma[i](x), x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := display(Graph, p)     end if   end do: end do:
 > # finite_element_xi is undefined finite_element_xi := Array(0..10, [seq(rand(0. .. 10.)(), k=1..11)]):   # it's up to you to provide the good values Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j));       p := plot(finite_element_sigma[i](x), x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := display(Graph, p)     end if   end do: end do:
 > # finite_element_epsilon is undefined finite_element_epsilon := Array(0..10, [seq(rand(0. .. 10.)(), k=1..11)]):   # it's up to you to provide the good values Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j));       p := plot(finite_element_sigma[i](x), x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := display(Graph, p)     end if   end do: end do:
 > # R(7)(x) is undefined R(7)(x)
 (1)
 > # You gave the expression of the R sunction, so why didn't use it? R := i -> E0*Heaviside(i+x0-x) + E1*Heaviside(i-(x+x0)): # example R(7); # But this likely won't give what you expect to get R(7)(x)
 (2)

First definition of  R

 > R := i -> E0*Heaviside(i+x0-x) + E1*Heaviside(i-(x+x0)): Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j));       p := plot(finite_element_sigma[i](x), x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := Graph, p     end if   end do: end do: plots:-display(Graph);
 > # do not write finite_element_sigma[i](x) Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j));       p := plot(finite_element_sigma[i], x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := Graph, p     end if   end do: end do: plots:-display(Graph);
 > # E10, E1, x0 are undefined E0 := 2:    # it's up to you to provide the good values E1 := 1:    # it's up to you to provide the good values x0 := 3:    # it's up to you to provide the good values Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j));       p := plot(finite_element_sigma[i], x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := Graph, p     end if   end do: end do: plots:-display(Graph);

Second definition of  R

 > R := i -> x -> E0*Heaviside(i+x0-x) + E1*Heaviside(i-(x+x0)): # example R(7)(x);
 (3)
 > Graph := NULL; for i from 0 to N-1 do   for j from 0 to 10 do     if j <= finite_element_xi[i] and finite_element_xi[i] <= j+1 then       finite_element_sigma[i] := evalf(finite_element_epsilon[i]*R(j)(x));       p := plot(finite_element_sigma[i], x = finite_element_xi[i] .. finite_element_xi[i+1]);       Graph := Graph, p     end if   end do: end do: plots:-display(Graph);
 >

## Explore...

Here is a Download Suggestion.mw based on the use of Explore Follow_me

A screen capture of what you get:

## Three ways...

Leading idea: find the zeroes of diff(s, phi):

 > restart; W:=f*mu*M**2*(1-sqrt(1+2*phi/(mu*M**2)))+(1-f)*(1+3*beta-(1+3*beta-3*beta*phi+beta*phi**2)*exp(phi))+g*nu*M**2*(1-sqrt(1-2*phi/(nu*M**2)))+(1-g)*(1+3*beta-(1+3*beta+3*beta*phi+beta*phi**2)*exp(-phi)): s:= subs(g=0.95,mu=1,beta=0.19,nu=0.5,M=2.277, f=0.4490, W): plot(s,phi=-2.3..0.6); # As it is seen for beta=0.19 and f=0.449 we have three consecutive local extrema. How can I find these critical values of (beta,f)?
 > ds    := diff(s, phi): omega := -2.3..0.6: plot(ds, phi=omega);
 > # Using solve: one solution is missing fnormal~({solve(ds)}); Locations := select((x -> verify(x, omega, 'interval')), %);
 (1)
 > # simplify(convert(ds, rational), size): # solve(%): # allvalues(%):
 > # Using fsolve: the method is governed by what you see on the plot above loc_1 := fsolve(ds, phi=omega); loc_2 := fsolve(ds, phi=loc_1..op(2, omega)); loc_3 := fsolve(ds, phi=op(1, omega)..loc_1); loc_4 := fsolve(ds, phi=op(1, omega)..loc_3-1e-6); Locations := fnormal~({loc_1, loc_2, loc_3, loc_4})
 (2)
 > # Using RootFinding:-NextZero: a safer and clever way than using fsolve fds := unapply(ds, phi): KeepSearching := true: start         := -2.3: RightBound    := 0.6: Locations     := NULL: while KeepSearching do   loc := RootFinding:-NextZero(phi -> fds(phi), start);   if loc < RightBound then     Locations     := Locations, loc:     start         := loc:   else     KeepSearching := false:   end if end do: Locations := fnormal~([Locations])
 (3)
 >