mmcdara

7891 Reputation

22 Badges

9 years, 54 days

MaplePrimes Activity


These are answers submitted by mmcdara

Probably a hideous way to proceed

Using  result := ssystem("tasklist <opts>"), where < opts >depends on what OS you use (Windows only) you can get a table (iin Windows sense) which contains all the PIDs.
This require to "postprocess" resultto display only the PIDs associated to a mserver.exeprocess

When you run xmaple, open a worksheet run the previous command to get the mserver's PIDs.
Proceed this way ach time you open a new worksheet: by comparing the new PID's list to the one of the previously opened worksheet you will get the PID associated to your current worksheet.
This can be useful is the stopbutton doesn't cancel the computations and if you can't close the current document.

Unfortunately I don't know how to identify to what node a mserverprocess belongs to when I distribute the computations over several nodes while using the GridPackage (being able to kill the computations on a specific node is something I would really like to have, all the more so that even after leaving Maple the process still exists and keeps using memory).

The command which generates output (4) is incorrect:

  • You use = instead of := 
  • You use exp^(...) instead of exp(...) 

Here is the correction (done with Maple 2015 whichh doesn't accept this awfull notation C(x, t) := ...)

 

restart

with(plottools):

with(plots):

with(CurveFitting):

Digits := 10:

f := proc (t) options operator, arrow; 7.0*exp(-(1/2000000)*(t-13180)^2)+4.7*exp(-(1/3200000)*(t-16000)^2) end proc:

p1 := plot(f(t), t = 0 .. 20000, color = green); -1; plots[display]({p1})

 

D1 := 15:

varepsilon := 200000:

L := 6500:

n := 200:

t := 1000

1000

(1)

lambda := simplify(evalf(n*Pi*sqrt((1/2)*D1+sqrt((1/4)*D1^2+varepsilon*(n*Pi/L)^2))/L))

.6928578233

(2)

b := 2*(int(f(t)*sin(m*Pi*x/L), x = 0 .. L))/L

-0.6366197724e-1*(0.1409730543e-28*cos(3.141592654*m)-0.1409730543e-28)/m

(3)

C := proc (x, t) options operator, arrow; sum(b*exp(-lambda^2*t)*sin(m*Pi*x/L), m = 1 .. 2) end proc

proc (x, t) options operator, arrow; sum(b*exp(-lambda^2*t)*sin(m*Pi*x/L), m = 1 .. 2) end proc

(4)

uu1000 := [seq(evalf(C(L-i, t)), i = 0 .. 6500, 100)]:

NULL

xx := [seq(k, k = 0 .. 6500, 100)]:

p2 := plot(xx, uu1000, color = cyan):

plots[plots:-display]({p2});

 
 

Download easy_way_mmcdara.mw

L:=[[3, 2], [2, 1], [1, 2], [1, 2], [2, 3], [2, 1], [1, 2], [1, 1], [2, 1], [1, 2], [1, 1], [2, 1], [1, 1], [1, 3], [1, 2], [2, 1], [1, 3], [1, 3], [1, 3], [1, 2], [2, 2], [2, 3]]: 

Ls := convert~(L0, set): 
L1 := map(u -> [[u[]], ListTools:-Occurrences(u, Ls)], {Ls[]}): 
Lk := map(u -> if numelems(u[1])=1 then [[u[1][1]$2], u[2]] else u end if, L1)

{[[1, 1], 3], [[1, 2], 11], [[1, 3], 4], [[2, 2], 1], [[2, 3], 3]}

L2 := [map(u -> [add(op(u[1])), u[2]], Lk)[]]

           [[2, 3], [3, 11], [4, 1], [4, 4], [5, 3]]

answer :=(`*`@op)( (`*`@op)~(L2))

                             190080


If you want to find the argument which minimizes the sum of the norms for a given value of lambda, this is quite simple:

w[opt] := (X^+ . X + lambda*~J)^(-1) . X^+ . Y

where X is a N by P matrix, Y a column vector of length N and J the identity matrix of dimension N.

But maybe you are talking about Ridge regression where you seek fot the couple (w[opt], lambda[opt]) which minimizes thz sum of the norms???

Several things are weird or incorrect:

  • Do you really need this weird Digits:=100???
    It works quite the same with Digits:=10.
     
  • p(t) is a nonsense, for the integration is done over t and thus its result cannot depend on t:
    p := int(f(t)*exp(-lambda^2*t), t=0..L)
                                         -134
                          -4.804438634 10    
    
  • As p is a constant, C(x, t) is a function of x alone.
    The message 
    Error, (in sum) summation variable previously assigned, second argument evaluates to 200. = 1 .. 500
    is clear: you have previously assigned n to 200.
    The cure: protect n by writing 'n' =1..500
    C := x -> sum((2/L*sin(n*Pi*x/L)*exp(1)*sin(n*Pi*x/L))*p, 'n' = 1 .. 500):
    C(1);
                                         -136
                          -1.871559740 10    
    

With these correction the worksheet doesn't contain any error...
But are you sure that the integral

int(f(t)*exp(-lambda^2*t), t=0..L)

is the one you wanted to write???



In Maple 2015.2 SumTools doesn't offer a procedure Change as IntegrationTools does.
So the idea to transform the sum/Sum into Int, to apply Change, and to go back to the sum/Sum.

S := sum(a[k]*(k+r)*(k+r-1)*x^(k+r-1), k = 0 .. infinity): 
J := Int(op(S)): 
IntegrationTools:-Change(J, k=j+1, j):
IntegrationTools:-Change(%, j=k, k): 
Sum(op(%))
         infinity                                     
          -----                                       
           \                                          
            )                                  (k + r)
           /     a[k + 1] (k + 1 + r) (k + r) x       
          -----                                       
          k = -1                                      

 

@sursumCorda 

You're right.
Here is another way based on the following principle:

  • Let A some term in rules.
  • Set B = eval(A, R) where R denotes the rules with element removed
  • Keep proceeding this way while B <> A after setting R=R \ B and A=B.

I initially removed elements of the form b=b (trivial 1-cycles).
 

restart

with(ListTools):

cycle := proc(i)
  local A, S, R, B, j:

  A := rules[i];
  S := A;
  if Search((rhs=lhs)(A), rules) > 0 then
    #----- cycle of length 2
    return [A, (rhs=lhs)(A)];
  end if:

  R := subsop(i=NULL, rules);
  B := lhs(A) = eval(rhs(A), R);

  if (lhs=rhs)(B) then
    return S, rhs(A)=rhs(B);

  else
    while numelems(R) > 0 and `not`(is(A=B)) do
     S := S, rhs(A)=rhs(B);
     j := Search(S[-1], R);
     A := R[j];
     R := subsop(j=NULL, R);

     if Search(rhs(A)=lhs(S[1]), rules) > 0 then
       #----- cycle of length 2
       return [S, rhs(A)=lhs(S[1])];
     end if:

     B := lhs(A) = eval(rhs(A), R);
    end do;
  end if:
  S := [S]:

  if lhs(S[1])=rhs(S[-1]) then
    if numelems(S) > 1 then
      return [S]
    else
      return []
    end if:
  else
      return []
  end if;
end proc:
 

rules := [g = d, e = a, a = b, f = d, c = g, b = b, d = c, c = f, a = e, f = g];
rules := map(u -> if `not`(is(u)) then u end if, rules);

[g = d, e = a, a = b, f = d, c = g, b = b, d = c, c = f, a = e, f = g]

 

[g = d, e = a, a = b, f = d, c = g, d = c, c = f, a = e, f = g]

(1)

for i from 1 to numelems(rules) do
  printf("from = %a : cycle = %a\n", rules[i], cycle(i));
end do;

from = g = d : cycle = [g = d, d = c, c = g]
from = e = a : cycle = [e = a, a = e]

from = a = b : cycle = []
from = f = d : cycle = [f = d, d = c, c = f]
from = c = g : cycle = [c = g, g = d, d = c]
from = d = c : cycle = [d = c, c = g, g = d]
from = c = f : cycle = [c = f, f = d, d = c]
from = a = e : cycle = [a = e, e = a]
from = f = g : cycle = [f = g, g = d, d = c, c = f]

 

 


 

Download cycles_1.mw


EDITED
to get all the cycles without permutations removed, replace the last for..do..end do structure by
 

# AEC stands for All Elementary Cycles (rotations removed) 

AEC := NULL:
for i from 1 to numelems(rules) do
  LocalCycle := cycle(i):
  if `not`(has(convert~([AEC], set), {convert(LocalCycle, set)})) then
    AEC := AEC, LocalCycle;   end if:   printf("from = %a : cycle = %a\n", rules[i], LocalCycle);
end do:
from = g = d : cycle = [g = d, d = c, c = g]
from = e = a : cycle = [e = a, a = e]
from = a = b : cycle = []
from = f = d : cycle = [f = d, d = c, c = f]
from = c = g : cycle = [c = g, g = d, d = c]
from = d = c : cycle = [d = c, c = g, g = d]
from = c = f : cycle = [c = f, f = d, d = c]
from = a = e : cycle = [a = e, e = a]
from = f = g : cycle = [f = g, g = d, d = c, c = f]

print~({AEC} minus {[]}):
                         [e = a, a = e]
                     [f = d, d = c, c = f]
                     [g = d, d = c, c = g]
                  [f = g, g = d, d = c, c = f]

 


The words you used made me think of your problem in terms of a graph.
The underlying idea is to represent your list of rules by a directed graph and then search for "directed cycles'.

I'm using Maple 2015.2 right now, which means I miss some recent features of GraphTheory (for instance single loops are nort allowed).
The code above finds:

  • 1-cycles, corresponding for instance to a rule such that a=a
  • 2-cycles, corresponding for instance to a couple of rules such that a=b, b=a
  • and FUNDAMENTAL cycles of higher lengths through the procedure CycleBasis (applied to the underlying undirected graph).

This means that some cycles may be missed, for instance those that could be created by a composition of cycles in the cycle basis.

Maybe someone will be interested in this alternative approach?

Download cycles.mw

@Carl Love is right.
But, one may infer the solution this way:

  • Replace 2014 by 2 (and ths 2015 by 3 and 2013 by 1)
  • Solving this 2 by 2 system gives x1=x2=1/2
  • First inference: could it be that x1=x2...=x2014=1/2014 ?
  • Second inference (used below): could it be that x1=x2...=x2014=x, wher "x" is to de determined
     

Answer:  x1 = x2 ... = x2014 = 1/2014

restart

sys := {
         sqrt(1+x)+sqrt(1+y)=2*sqrt(3/2),
         sqrt(1-x)+sqrt(1-y)=2*sqrt(1/2)
       };

{(1-x)^(1/2)+(1-y)^(1/2) = 2^(1/2), (1+x)^(1/2)+(1+y)^(1/2) = 6^(1/2)}

(1)

solve(sys)

{x = 1/2, y = 1/2}

(2)

solve(2014*sqrt(1+x)=2014*sqrt(2015/2014));
solve(2014*sqrt(1-x)=2014*sqrt(2013/2014));

1/2014

 

1/2014

(3)

 

Download InferedSolution.mw

The x label is conform to what you want but the construction is quite complex (it uses MathML).

  • "&#x28;" codes the left parenthesis
  • "&#x29;" codes the right one
  •  mo("-",mathcolor="white") is a trick to introduce a space
  •  msup(mo("x"),mo("n")) generates xn

Don't focus ob what the y label looks like: you need to execute the code to get a correct display.

restart:

with(Units):

plot(
  x, x=0..1,
  labels=[
    typeset(`#mi("A")`, `#mrow(mo("-",mathcolor="white"),mo("&#x28;"),msup(mo("K/m"),mo("3")),mo("&#x29;"))`),
    A*Unit('K'/'m^3')
  ]
)

 

 

Download TwoWays.mw

For your system to be a PDE system it has to contain partial derivatives wrt both x AND t.
As no derivative wrt t is present, you have a simple ODE system.

You already noticed that writting eq3 and eq4, and then did right solving this system.
If what puzzles you is that the solution doesn't depend on arbtrary functions of t you can do somethinh like that:
in the file below f(t) and g(t) are arbitrary functions of t, those who Maple would name _F1(t) and _F2(t)
 

restart

eq1 := diff(p(x, t), x) = p(x, t)*(1/(p(x, t)^2*q(x, t)^(1/3)))

diff(p(x, t), x) = 1/(p(x, t)*q(x, t)^(1/3))

(1)

eq2 := diff(q(x, t), x) = -(-3*q(x, t))*(1/(p(x, t)^2*q(x, t)^(1/3)))

diff(q(x, t), x) = 3*q(x, t)^(2/3)/p(x, t)^2

(2)

sys := {eq1, eq2}:

SYS := eval(sys, {p(x, t)=P(x), q(x, t)=Q(x)});

{diff(P(x), x) = 1/(P(x)*Q(x)^(1/3)), diff(Q(x), x) = 3*Q(x)^(2/3)/P(x)^2}

(3)

SOL := dsolve(SYS);

[{Q(x) = _C1*x+_C2}, {P(x) = 3^(1/2)*((diff(Q(x), x))*Q(x)^(2/3))^(1/2)/(diff(Q(x), x)), P(x) = -3^(1/2)*((diff(Q(x), x))*Q(x)^(2/3))^(1/2)/(diff(Q(x), x))}]

(4)

sol := eval(SOL, {P(x)=p(x, t), Q(x)=q(x, t)});

[{q(x, t) = _C1*x+_C2}, {p(x, t) = 3^(1/2)*((diff(q(x, t), x))*q(x, t)^(2/3))^(1/2)/(diff(q(x, t), x)), p(x, t) = -3^(1/2)*((diff(q(x, t), x))*q(x, t)^(2/3))^(1/2)/(diff(q(x, t), x))}]

(5)

sol := eval(sol, {_C1=f(t), _C2=g(t)}):

# check if the initial system is verified for the first solution

sol_1 := [op(sol[1]), eval(op(sol[2])[1], op(sol[1]))];
eval([sys[]], sol_1);
is~(%)

[q(x, t) = f(t)*x+g(t), p(x, t) = 3^(1/2)*(f(t)*(f(t)*x+g(t))^(2/3))^(1/2)/f(t)]

 

[(1/3)*f(t)*3^(1/2)/((f(t)*(f(t)*x+g(t))^(2/3))^(1/2)*(f(t)*x+g(t))^(1/3)) = (1/3)*f(t)*3^(1/2)/((f(t)*(f(t)*x+g(t))^(2/3))^(1/2)*(f(t)*x+g(t))^(1/3)), f(t) = f(t)]

 

[true, true]

(6)

# check if the initial system is verified for the second solution

sol_2 := [op(sol[1]), eval(op(sol[2])[2], op(sol[1]))];
eval([sys[]], sol_1);
is~(%)

[q(x, t) = f(t)*x+g(t), p(x, t) = -3^(1/2)*(f(t)*(f(t)*x+g(t))^(2/3))^(1/2)/f(t)]

 

[-(1/3)*f(t)*3^(1/2)/((f(t)*(f(t)*x+g(t))^(2/3))^(1/2)*(f(t)*x+g(t))^(1/3)) = -(1/3)*f(t)*3^(1/2)/((f(t)*(f(t)*x+g(t))^(2/3))^(1/2)*(f(t)*x+g(t))^(1/3)), f(t) = f(t)]

 

[true, true]

(7)

# complete solution

print~(sol_1):
print(cat("-"$80)):
print~(sol_2):

q(x, t) = f(t)*x+g(t)

 

p(x, t) = -3^(1/2)*(f(t)*(f(t)*x+g(t))^(2/3))^(1/2)/f(t)

 

"--------------------------------------------------------------------------------"

 

q(x, t) = f(t)*x+g(t)

 

p(x, t) = -3^(1/2)*(f(t)*(f(t)*x+g(t))^(2/3))^(1/2)/f(t)

(8)

 

Download Solve_as_an_ODE_system.mw

When I get some "counter-contuitive" result, I use to try and figure out what could have happened, for instance by seeking to the problem otherwise (please don't take this as a critic).
In your case I would have done this:

  • Seek if there exist real solutions by checking if the operands under each square root are positive on a non empty interval.
  • If not, replace x by y*I*z and operate explicitely in the  C field.
  • Observe the solutions.
    Evaluate y*I*z for these solutions.

You will see that the result of this last step (1/9 and 2) is exactely the result that Maple gave you.
So I guess thet Maple takes the same path that I described above (?)

restart

a1 := sqrt(x^2-10*x+1):
a2 := sqrt(-8*x^2+9*x-1):

s1 := solve(op(1, a1) >= 0);
s2 := solve(op(1, a2) >= 0);

RealRange(-infinity, 5-2*6^(1/2)), RealRange(5+2*6^(1/2), infinity)

 

RealRange(1/8, 1)

(1)

infolevel[solve] := 4:
s12 := solve({op(1, a1) >= 0, op(1, a2) >= 0} );

Main: Entering solver with 2 equations in 1 variable
Main: attempting to solve as a linear system
Main: system cannot be directly solved as a linear system
Main: attempting to solve as a semi-algebraic system

Main: Semi-algebraic solver successful. Exiting solver returning 1 solution
solve: Warning: no solutions found

 

# There is thus no real root to f


f:= a1-a2:
plots:-complexplot(f, x=-1..1, gridlines=true, numpoints=400)

 

# as there exist no real value of sucg tht f=0,
# let's write x=y+z*I

If := eval(f, x = y+I*z)

((y+I*z)^2-10*y-(10*I)*z+1)^(1/2)-(-8*(y+I*z)^2+9*y+(9*I)*z-1)^(1/2)

(2)

sol := solve(If)

Main: Entering solver with 1 equation in 2 variables

Main: attempting to solve as a linear system

Dispatch: dispatching to Radicals handler

Transformer:   solving for linear equation in _S000001
Recurse: recursively solving 1 equations in 2 variables
Dispatch: handling a single polynomial
Main: polynomial system split into 1 parts under preprocessing
Main: subsystem is essentially univariate
UnivariateHandler: subsystem has only one equation

UnivariateHandler: subsystem has only one equation
UnivariateHandler: solving as if univariate in y
Polynomial: solving a linear polynomial
UnivariateHandler: subsystem has only one equation
UnivariateHandler: solving as if univariate in y
Polynomial: solving a linear polynomial
Main: solving successful - now forming solutions
Main: Exiting solver returning 2 solutions

 

{y = -I*z+2, z = z}, {y = -I*z+1/9, z = z}

(3)

eval(If, sol[1]);
eval(If, sol[2]);

0

 

0

(4)

 

Download Another_POV.mw

Finally:

eval(y+I*z, sol[1]), eval(y+I*z, sol[2]);
                                 1
                              2, -
                                 9

 

If you really want the term within the parenthesis to  be of the form xn-1*((n+1)*x-2n), you can do this

map(collect, simplify(eval(f, x^n='x'*x^(n-1)), size), x)

                                      (n - 1)
                   ((n + 1) x - 2 n) x       

"Proof" is within the worksheet

 

restart

a1 := 1;

1

(1)

g[2] := -S0*eta[2]+b+r;

-S0*eta[2]+b+r

(2)

NULL

a[1] := c+epsilon+g[2];

-S0*eta[2]+b+c+epsilon+r

(3)

a[2] := epsilon*g[2]+c*(epsilon+g[2])-S0*b*eta[3];

epsilon*(-S0*eta[2]+b+r)+c*(-S0*eta[2]+b+epsilon+r)-S0*b*eta[3]

(4)

a[3] := -S0*Pi*b*eta[1]+c*(-S0*b*eta[3]+epsilon*g[2])

-S0*Pi*b*eta[1]+c*(-S0*b*eta[3]+epsilon*(-S0*eta[2]+b+r))

(5)

eq := lambda^3+lambda^2*a[1]+lambda*a[2]+a[3]

lambda^3+lambda^2*(-S0*eta[2]+b+c+epsilon+r)+lambda*(epsilon*(-S0*eta[2]+b+r)+c*(-S0*eta[2]+b+epsilon+r)-S0*b*eta[3])-S0*Pi*b*eta[1]+c*(-S0*b*eta[3]+epsilon*(-S0*eta[2]+b+r))

(6)

# eq has either 1 or 3 real roots.

with(LargeExpressions):
collect(eval(eq, lambda=u-a[1]/3), u, Veil[C]);

R := 4* (-1/3*C[1]) + 27 * (-1/27*C[2])^2:

Unveil[C](R):

R := simplify(%, size);




 

u^3-(1/3)*C[1]*u-(1/27)*C[2]

 

-(4/3)*S0^2*eta[2]^2+(8/3)*S0*b*eta[2]-4*S0*b*eta[3]-(4/3)*S0*c*eta[2]-(4/3)*S0*epsilon*eta[2]+(8/3)*S0*r*eta[2]-(4/3)*b^2+(4/3)*b*c+(4/3)*b*epsilon-(8/3)*b*r-(4/3)*c^2+(4/3)*c*epsilon+(4/3)*c*r-(4/3)*epsilon^2+(4/3)*epsilon*r-(4/3)*r^2+(4/27)*(b^3+((-3*eta[2]+(9/2)*eta[3])*S0-(3/2)*c+3*r-(3/2)*epsilon)*b^2+(3*eta[2]*(eta[2]-(3/2)*eta[3])*S0^2+((3*eta[2]-9*eta[3])*c+(-6*eta[2]+(9/2)*eta[3])*r+(3*eta[2]+(9/2)*eta[3])*epsilon-(27/2)*Pi*eta[1])*S0-(3/2)*c^2+(-3*r+6*epsilon)*c-3*epsilon*r+3*r^2-(3/2)*epsilon^2)*b+((1/2)*S0*eta[2]+c-(1/2)*r-(1/2)*epsilon)*(-S0*eta[2]+c-2*epsilon+r)*(2*S0*eta[2]+c+epsilon-2*r))^2

(7)

# R > 0 ==> only one real root
# R = 0 ==> a double real root and a simple real root
# R < 0 ==> three different real roots
#
# R depends on 8 parameters

params := indets(R, name) minus {Pi};

# And so there ar all the chances at world for R to be a non positive or a non negative
# function, which means that one can probably find stes of values for these 8 parameters
# such that R <, or R=0, or R > 0.
# For instance

collect(R, eta[1], Veil[K]);  # R being a degree 2 polynomial wrt eta[1] there exists potentially
                              # two real values of eta[1] such that R=0 for any choice of the 7
                              # remaining parameters.
                              # Let u the lower root and v the upper one (assuming they are both
                              # real: depending on the sign of K[1], R > 0 (or R < 0) within
                              # (u, v) and R < 0 (or R > 0) within (-oo, u) union (v, +oo).
                              #
                              # Which means that having set to numerical values all the parameters
                              # but eta[1] can potentially lead to R < 0, or R = 0, or R > 0...
                              # ... unless v is always less than 0 in which case, as eta[1] is
                              # assumed positive, R has always the same sign.
 

{S0, b, c, epsilon, r, eta[1], eta[2], eta[3]}

 

27*K[1]*eta[1]^2+2*K[2]*eta[1]+(1/27)*K[3]

(8)

# the analysis is far more complex several other parameters for R is then a 6th degree polynomial
# equation of any of them.

print~( map(u -> [u, collect(R, u, Veil[K])], params) ):

[S0, (4/27)*K[11]*S0^6-(4/9)*K[12]*S0^5+(1/9)*K[13]*S0^4-(2/27)*K[14]*S0^3+(1/9)*K[15]*S0^2-(2/9)*K[16]*S0+(1/27)*K[17]]

 

[b, (4/27)*b^6-(4/9)*K[18]*b^5+(1/9)*K[19]*b^4-(2/27)*K[20]*b^3+(1/9)*K[21]*b^2-(2/9)*K[22]*b+(1/27)*K[23]]

 

[c, (4/27)*c^6+(4/9)*K[24]*c^5-(1/9)*K[25]*c^4-(2/27)*K[26]*c^3-(1/9)*K[27]*c^2+(2/9)*K[28]*c+(1/27)*K[29]]

 

[epsilon, (4/27)*epsilon^6+(4/9)*K[30]*epsilon^5-(1/9)*K[31]*epsilon^4-(2/27)*K[32]*epsilon^3-(1/9)*K[33]*epsilon^2+(2/9)*K[34]*epsilon+(1/27)*K[35]]

 

[r, (4/27)*r^6-(4/9)*K[36]*r^5+(1/9)*K[37]*r^4-(2/27)*K[38]*r^3+(1/9)*K[39]*r^2-(2/9)*K[40]*r+(1/27)*K[41]]

 

[eta[1], 27*K[1]*eta[1]^2+2*K[2]*eta[1]+(1/27)*K[3]]

 

[eta[2], (4/27)*K[4]*eta[2]^6-(4/9)*K[5]*eta[2]^5+(1/9)*K[6]*eta[2]^4-(2/27)*K[7]*eta[2]^3+(1/9)*K[8]*eta[2]^2-(2/9)*K[9]*eta[2]+(1/27)*K[10]]

 

[eta[3], 3*K[42]*eta[3]^2+(2/3)*K[43]*eta[3]+(1/27)*K[44]]

(9)

# then is completely impossible to assess the number of real lambda roots not even
# saying their signs

Download roots.mw

If I understand correctly, here is a solution
(I sorted the selected columns: if you don't want to sort them modify my code accordingly)

restart:

# Let A some N by M matrix, for instance

A := (N, M) -> Matrix(N, M, (i, j) -> i+N*j):

ChooseColumns := proc(Mat::Matrix, k::posint)
  local m, n, choice:
  n, m := LinearAlgebra:-Dimensions(Mat):
  if k > m-1 then
    error  cat("k (", k, ") must be lower the the number of tows (", m, ")")
  end if:
  choice := sort(combinat:-randperm([$1..m-1])[1..k]):
  return Mat[.., [choice[], m]]
end proc:

B := A(10, 5):
B, ChooseColumns(B, 3);

Matrix(10, 5, {(1, 1) = 11, (1, 2) = 21, (1, 3) = 31, (1, 4) = 41, (1, 5) = 51, (2, 1) = 12, (2, 2) = 22, (2, 3) = 32, (2, 4) = 42, (2, 5) = 52, (3, 1) = 13, (3, 2) = 23, (3, 3) = 33, (3, 4) = 43, (3, 5) = 53, (4, 1) = 14, (4, 2) = 24, (4, 3) = 34, (4, 4) = 44, (4, 5) = 54, (5, 1) = 15, (5, 2) = 25, (5, 3) = 35, (5, 4) = 45, (5, 5) = 55, (6, 1) = 16, (6, 2) = 26, (6, 3) = 36, (6, 4) = 46, (6, 5) = 56, (7, 1) = 17, (7, 2) = 27, (7, 3) = 37, (7, 4) = 47, (7, 5) = 57, (8, 1) = 18, (8, 2) = 28, (8, 3) = 38, (8, 4) = 48, (8, 5) = 58, (9, 1) = 19, (9, 2) = 29, (9, 3) = 39, (9, 4) = 49, (9, 5) = 59, (10, 1) = 20, (10, 2) = 30, (10, 3) = 40, (10, 4) = 50, (10, 5) = 60}), Matrix(10, 4, {(1, 1) = 11, (1, 2) = 21, (1, 3) = 31, (1, 4) = 51, (2, 1) = 12, (2, 2) = 22, (2, 3) = 32, (2, 4) = 52, (3, 1) = 13, (3, 2) = 23, (3, 3) = 33, (3, 4) = 53, (4, 1) = 14, (4, 2) = 24, (4, 3) = 34, (4, 4) = 54, (5, 1) = 15, (5, 2) = 25, (5, 3) = 35, (5, 4) = 55, (6, 1) = 16, (6, 2) = 26, (6, 3) = 36, (6, 4) = 56, (7, 1) = 17, (7, 2) = 27, (7, 3) = 37, (7, 4) = 57, (8, 1) = 18, (8, 2) = 28, (8, 3) = 38, (8, 4) = 58, (9, 1) = 19, (9, 2) = 29, (9, 3) = 39, (9, 4) = 59, (10, 1) = 20, (10, 2) = 30, (10, 3) = 40, (10, 4) = 60})

(1)

 

Download ChooseColumns.mw

BTW, did you succeed in loading the PLS_PCR_Linnerund.mw file ? 
(see my las t replyy to the related thread).

First 26 27 28 29 30 31 32 Last Page 28 of 65