dharr

Dr. David Harrington

8220 Reputation

22 Badges

20 years, 339 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

Indexed notation like cov[x]:=newval; would be the usual way to modify an entry. I don't see anything wrong with your code; I'd use eval(cov) over op(cov) for the return value for a slight improvement in readability. To reset a table, just set it to a new blank table: cov:=table();

In Maple 2023, cov := table([seq(x = x, x = X)]) could be cov := table(X =~ X)

and cov := table([seq(x = NULL, x = X)]) could be cov:=table(X =~ NULL)

Negative values (as opposed to index=negative value) are just estimates of roots. index and allvalues work best for polynomials. I do it also by fsolve, by you can also use NextZero for systematically finding roots.

I don't know another way to add "index=number", other than convert(r1,RootOf,form=index) which gives index=1 (for a polynomial).

restart

arg := _Z*cos(_Z)-sqrt(sin(_Z)^2)

_Z*cos(_Z)-(sin(_Z)^2)^(1/2)

r0 := RootOf(arg)

RootOf(_Z*cos(_Z)-(sin(_Z)^2)^(1/2))

inter := -6 .. 6

-6 .. 6

plot(arg, _Z = inter, size = [300, 300])

Get any one root

evalf(r0); fsolve(arg); fsolve(arg, _Z = inter)

0.

0.

-2.028757838

If you are willing to look at the plot and do hand estimates then RootOf (or fsolve) works. Save this input line first, and then evaluate.

vals := 'evalf(RootOf(arg, -4)), evalf(RootOf(arg, -2)), evalf(RootOf(arg, 0)), evalf(RootOf(arg, 5))'

evalf(RootOf(arg, -4)), evalf(RootOf(arg, -2)), evalf(RootOf(arg, 0)), evalf(RootOf(arg, 5))

vals

-4.493409458, -2.028757838, 0., 4.913180439

Or intervals are probably more reliable.

evalf(RootOf(arg, -5 .. -4)), evalf(RootOf(arg, -3 .. -1)), evalf(RootOf(arg, -1 .. 1)), evalf(RootOf(arg, 4 .. 5))

-4.493409458, -2.028757838, 0., 4.913180439

Conversions are possible - the strange result here is because index is really  only for polynomials

r1 := RootOf(arg, -4); convert(r1, RootOf, form = index); evalf(%)

RootOf(_Z*cos(_Z)-(sin(_Z)^2)^(1/2), -4)

RootOf(_Z*cos(_Z)-(eval(RootOf(_Z^2-E, index = 1), {E = sin(_Z)^2})), -4)

-4.493409458

Not sure whats happening in the next two

convert(r1, 'RootOf', form = 'interval')

Error, (in type/algext) too many levels of recursion

convert(RootOf(arg, -5 .. -4), RootOf, form = numeric)

Error, (in type/algext) too many levels of recursion

Using fsolve and avoid

a1 := fsolve(arg, _Z = inter); a2 := fsolve(arg, _Z, avoid = {_Z = a1}); a3 := fsolve(arg, _Z, avoid = {_Z = a1, _Z = a2}); a4 := fsolve(arg, _Z, avoid = {_Z = a1, _Z = a2, _Z = a3})

-2.028757838

0.

-4.493409458

4.913180439

NULL

Download RootOf.mw

In future it would be helpful if you uploaded your worksheet with the green up-arrow in the Mapleprimes editor. One problem is that your code shows diff*(S(t), t) not diff(S(t), t) because you have a space between diff and (, which is interpreted as multiplication. (also in other places).

Here is a solution, based on iterating through the related binary trees - see https://en.wikipedia.org/wiki/Catalan_number

restart

Consider a regular polygon with m = n+2 vertices. We want all non-crossing triangulations of the vertices, see https://en.wikipedia.org/wiki/Catalan_number 
These are in 1:1 correspondence with the binary trees with n+1 leaves or n internal nodes, or with n deep nested parentheses.
Consider the example of a pentagon.

with(GraphTheory)

m := 5; n := m-2

5

3

The edges of the pentagon are labelled A,B,C,D,1. 1 corresponds to the first level vertex of the binary tree, and the others to the leaves. Edges 2..n will be internal edges added, which correspond to the internal vertices of the binary tree.

leaves := [seq({i, i+1}, i = 1 .. m-1)]; edges := {leaves[], {1, m}}; G := Graph(edges)

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

{{1, 2}, {1, 5}, {2, 3}, {3, 4}, {4, 5}}

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 3}, (3) = {2, 4}, (4) = {3, 5}, (5) = {1, 4}}), `GRAPHLN/table/1`, 0)

The P representation (see ?Iterator:-Trees) for binary trees/nested parentheses gives the internal tree vertices, which are interspersed with the leafs to give a nested parentheses representation,

e.g. [1,3,2] means A 1 B 3 C 2 D or  A1((B3C)2D) or just A((BC)D), where the 3 indicates the first pairing and 2 the second and 1 the last.
Starting with 3, we generate internal edge 3 as the third side of the triangle B3C, then internal edge 2 is the third side of the triangle 32D, then the last triangle is 21A

p := Array([1, 3, 2])

Array(%id = 36893489648404452884)

The following (inelegant and inefficient) code triangulates this case, contracting the parentheses in turn

plist := convert(p, list);
edgelist := leaves;
for i from n by -1 to 2 do
    member(i, plist, 'j');
    print(i, j);
    newedge := symmdiff(edgelist[j], edgelist[j + 1]);
    AddEdge(G, newedge);
    edgelist := [edgelist[1 .. j - 1][], newedge, edgelist[j + 2 .. -1][]];
    plist := subsop(j = NULL, plist);
end do;

[1, 3, 2]

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

true

3, 2

{2, 4}

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 3, 4}, (3) = {2, 4}, (4) = {2, 3, 5}, (5) = {1, 4}}), `GRAPHLN/table/1`, 0)

[{1, 2}, {2, 4}, {4, 5}]

[1, 2]

true

2, 2

{2, 5}

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 3, 4, 5}, (3) = {2, 4}, (4) = {2, 3, 5}, (5) = {1, 2, 4}}), `GRAPHLN/table/1`, 0)

[{1, 2}, {2, 5}]

[1]

DrawGraph(G, size = [200, 200])

So here is a procedure that iterates though all cases

triangulations:=proc(m::posint)
  local n,leaves,edges,i,t,j,G,T,p,plist,newedge,graphs,edgelist,ng;
  uses GraphTheory;
  n:=m-2;
  leaves := [seq({i, i + 1}, i = 1 .. m - 1)];
  edges := {leaves[], {1, m}};
  T:=Iterator:-BinaryTrees(n);
  graphs:=table();
  ng:=0;
  for t in T do
    G := Graph(convert(edges, set));
    p := Iterator:-Trees:-Convert(t, 'LR', 'P');
    plist := convert(p, list);
    edgelist := leaves;
    for i from n by -1 to 2 do
      member(i, plist, 'j');
      newedge := symmdiff(edgelist[j], edgelist[j + 1]);
      AddEdge(G, newedge);
      edgelist := [edgelist[1 .. j - 1][], newedge, edgelist[j + 2 .. -1][]];
      plist := subsop(j = NULL, plist);
    end do;
    ng:=ng+1;
    graphs[ng]:=copy(G);
  end do;
  convert(graphs,list);
end proc:

And here are the 14 cases for m=6

DrawGraph(triangulations(6))

 

 

 

 

 

NULL

Download Catalan.mw

Does work in 2-D

sequenceoperation.mw

Well that paper (published here) just deals with many special cases, which are straightforward, but tedious, to program in Maple. For example, here is the algorithm for quintics (I didn't deal with the cases with linear factors).

restart

Algorithm 4.1 from https://facstaff.elon.edu/cawtrey/acj-reducible.pdf

quinticGG:=proc(p,x::name)
local f,degrees,f1,f2,d1,d2;
if not type(p,polynom(integer,x)) then error "not a polynomial in %1 with integer coefficients",x end if;
if degree(p,x)<> 5 then error "polynomial not a quintic" end if;
if irreduc(p) then return GroupTheory:-GaloisGroup(p,x) end if;
f:=factors(p)[2];
degrees:=map(q->degree(q[1],x),f);
if has(degrees,1) then error "polynomial has linear factors" end if;
if degrees<>[2,3] and degrees<>[3,2] then error "degrees %1 unexpected",degrees end if;
if degrees=[2,3] then f1,f2:=map(q->q[1],f)[] else f2,f1:=map(q->q[1],f)[] end if;
d1:=discrim(f1,x);
d1:=d1*signum(d1);
d2:=discrim(f2,x);
d2:=d2*signum(d2);
if type(sqrt(d2),integer) then
  return GroupTheory:-CyclicGroup(6)
elif type(sqrt(d1*d2),integer) then
  return GroupTheory:-SymmetricGroup(3)
else
  return GroupTheory:-DihedralGroup(6)
end if;
end proc:

quinticGG(x^5-x+15, x)

_m1708329273568

quinticGG((x^2+3)*(x^3+2), x)

_m1708327735136

quinticGG((x^2+1)*(x^3-3*x+1), x)

_m1708326740160

quinticGG((x^2+1)*(x^3+2), x)

_m1708324937024

quinticGG(x^5-x+1, x)

_m1708323778880

NULL

Download GaloisAlgorithm.mw

ListTools:-Categorize is the easiest way to do this.

restart

with(GraphTheory)

graphs := [seq(RandomGraphs:-RandomGraph(5, 4), 1 .. 10)]

indextable := table(`~`[`=`](graphs, [`$`(1 .. nops(graphs))]))

classes := [ListTools:-Categorize(IsIsomorphic, graphs)]

[[GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {3, 5}, (2) = {}, (3) = {1, 5}, (4) = {5}, (5) = {1, 3, 4}}), `GRAPHLN/table/1`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 5}, (3) = {5}, (4) = {}, (5) = {1, 2, 3}}), `GRAPHLN/table/2`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 5}, (2) = {1, 4, 5}, (3) = {}, (4) = {2}, (5) = {1, 2}}), `GRAPHLN/table/4`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 4, 5}, (2) = {1, 5}, (3) = {}, (4) = {1}, (5) = {1, 2}}), `GRAPHLN/table/5`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {4, 5}, (2) = {5}, (3) = {}, (4) = {1, 5}, (5) = {1, 2, 4}}), `GRAPHLN/table/7`, 0)], [GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {4, 5}, (2) = {3}, (3) = {2}, (4) = {1, 5}, (5) = {1, 4}}), `GRAPHLN/table/3`, 0)], [GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 4}, (2) = {1, 3}, (3) = {2, 4}, (4) = {1, 3}, (5) = {}}), `GRAPHLN/table/6`, 0)], [GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 4}, (2) = {1}, (3) = {5}, (4) = {1, 5}, (5) = {3, 4}}), `GRAPHLN/table/8`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {4, 5}, (2) = {4}, (3) = {5}, (4) = {1, 2}, (5) = {1, 3}}), `GRAPHLN/table/9`, 0)], [GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {4}, (2) = {5}, (3) = {5}, (4) = {1, 5}, (5) = {2, 3, 4}}), `GRAPHLN/table/10`, 0)]]

uniquegraphs := map2(op, 1, classes)

[GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {3, 5}, (2) = {}, (3) = {1, 5}, (4) = {5}, (5) = {1, 3, 4}}), `GRAPHLN/table/1`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {4, 5}, (2) = {3}, (3) = {2}, (4) = {1, 5}, (5) = {1, 4}}), `GRAPHLN/table/3`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 4}, (2) = {1, 3}, (3) = {2, 4}, (4) = {1, 3}, (5) = {}}), `GRAPHLN/table/6`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {2, 4}, (2) = {1}, (3) = {5}, (4) = {1, 5}, (5) = {3, 4}}), `GRAPHLN/table/8`, 0), GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(1..5, {(1) = {4}, (2) = {5}, (3) = {5}, (4) = {1, 5}, (5) = {2, 3, 4}}), `GRAPHLN/table/10`, 0)]

map(proc (i) options operator, arrow; indextable[i] end proc, uniquegraphs)

[1, 3, 6, 8, 10]

seq(map(proc (i) options operator, arrow; indextable[i] end proc, class), `in`(class, classes))

[1, 2, 4, 5, 7], [3], [6], [8, 9], [10]

NULL

Download Categorize.mw

Since text files are sequential and not random access, I think you have to live with this inefficiency, and just read and discard the first lines, say

filnam:=cat(currentdir(),"/graph8c.txt");
to 199 while readline(filnam) <> 0 do end do:
for i to 101 while (gr[i]:=readline(filnam))<>0 do end do:

and then ConvertGraph to get them in Maple Graph form. (The end of file conditions should probably be better handled here.)

You could read the whole file in at once with readdata or FileTools:-Text:-ReadFile if memory isn't a limitation.

I'd use map here:

map(x->`if`(x<0.2,0,x),a);

 

The index=real[2] is the problem. If this is changed to a regular index (index=1 in this case), then Minpoly works as expected.

restart;

alias(`~`[`=`](alpha__ || (1 .. 3), ` $`, RootOf(11*_Z^9+17*_Z^8-64*_Z^7-280*_Z^6+142*_Z^5+370*_Z^4+376*_Z^3-96*_Z^2+47*_Z-11, .2246 .. .2266), RootOf(11*_Z^9+17*_Z^8-64*_Z^7-280*_Z^6+142*_Z^5+370*_Z^4+376*_Z^3-96*_Z^2+47*_Z-11, 1.671 .. 1.68), RootOf(11*_Z^9+17*_Z^8-64*_Z^7-280*_Z^6+142*_Z^5+370*_Z^4+376*_Z^3-96*_Z^2+47*_Z-11, 2.648 .. 2.657)))

alpha__1, alpha__2, alpha__3

(1)

(Not all of these are independent; we could have chosen just one.)

evala(Algfield({`&alpha;__1`, `&alpha;__2`, `&alpha;__3`}))

({PDETools:-Solve})({`~`[`>=`](a, b, ` $`, 0), a^5*b+4*a^4*b^2+4*a^3*b^3-7*a^4*b-6*a^2*b^3-7*a*b^4+b^5-6*a^3*b+12*a^2*b^2+4*b^4+4*a^3-6*a*b^2+4*b^3+4*a^2-7*a*b+a = 0, a <> b})

{{a = alpha__1, b = RootOf(1216*_Z^4+(264*alpha__1^8+408*alpha__1^7-1580*alpha__1^6-6832*alpha__1^5+3508*alpha__1^4+9944*alpha__1^3+9948*alpha__1^2-10752*alpha__1+5204)*_Z^3+(891*alpha__1^8+1652*alpha__1^7-4748*alpha__1^6-24076*alpha__1^5+5354*alpha__1^4+35356*alpha__1^3+29668*alpha__1^2-196*alpha__1+3971)*_Z^2+(506*alpha__1^8+980*alpha__1^7-2264*alpha__1^6-12420*alpha__1^5+3676*alpha__1^4+11596*alpha__1^3+33800*alpha__1^2-7772*alpha__1+1210)*_Z-473*alpha__1^8-720*alpha__1^7+2560*alpha__1^6+10960*alpha__1^5-8034*alpha__1^4-13840*alpha__1^3-9304*alpha__1^2+1104*alpha__1-1133, index = real[2])}, {a = alpha__2, b = RootOf(1216*_Z^4+(264*alpha__2^8+408*alpha__2^7-1580*alpha__2^6-6832*alpha__2^5+3508*alpha__2^4+9944*alpha__2^3+9948*alpha__2^2-10752*alpha__2+5204)*_Z^3+(891*alpha__2^8+1652*alpha__2^7-4748*alpha__2^6-24076*alpha__2^5+5354*alpha__2^4+35356*alpha__2^3+29668*alpha__2^2-196*alpha__2+3971)*_Z^2+(506*alpha__2^8+980*alpha__2^7-2264*alpha__2^6-12420*alpha__2^5+3676*alpha__2^4+11596*alpha__2^3+33800*alpha__2^2-7772*alpha__2+1210)*_Z-473*alpha__2^8-720*alpha__2^7+2560*alpha__2^6+10960*alpha__2^5-8034*alpha__2^4-13840*alpha__2^3-9304*alpha__2^2+1104*alpha__2-1133, index = real[2])}, {a = alpha__3, b = RootOf(1216*_Z^4+(264*alpha__3^8+408*alpha__3^7-1580*alpha__3^6-6832*alpha__3^5+3508*alpha__3^4+9944*alpha__3^3+9948*alpha__3^2-10752*alpha__3+5204)*_Z^3+(891*alpha__3^8+1652*alpha__3^7-4748*alpha__3^6-24076*alpha__3^5+5354*alpha__3^4+35356*alpha__3^3+29668*alpha__3^2-196*alpha__3+3971)*_Z^2+(506*alpha__3^8+980*alpha__3^7-2264*alpha__3^6-12420*alpha__3^5+3676*alpha__3^4+11596*alpha__3^3+33800*alpha__3^2-7772*alpha__3+1210)*_Z-473*alpha__3^8-720*alpha__3^7+2560*alpha__3^6+10960*alpha__3^5-8034*alpha__3^4-13840*alpha__3^3-9304*alpha__3^2+1104*alpha__3-1133, index = real[2])}}

(2)

bSol := `~`[subs](%, b)

evalf[2*Digits](`~`[eval](11*_X^9-47*_X^8+96*_X^7-376*_X^6-370*_X^5-142*_X^4+280*_X^3+64*_X^2-17*_X-11, `~`[`=`](_X, bSol)))

{RootOf(1216*_Z^4+(264*alpha__1^8+408*alpha__1^7-1580*alpha__1^6-6832*alpha__1^5+3508*alpha__1^4+9944*alpha__1^3+9948*alpha__1^2-10752*alpha__1+5204)*_Z^3+(891*alpha__1^8+1652*alpha__1^7-4748*alpha__1^6-24076*alpha__1^5+5354*alpha__1^4+35356*alpha__1^3+29668*alpha__1^2-196*alpha__1+3971)*_Z^2+(506*alpha__1^8+980*alpha__1^7-2264*alpha__1^6-12420*alpha__1^5+3676*alpha__1^4+11596*alpha__1^3+33800*alpha__1^2-7772*alpha__1+1210)*_Z-473*alpha__1^8-720*alpha__1^7+2560*alpha__1^6+10960*alpha__1^5-8034*alpha__1^4-13840*alpha__1^3-9304*alpha__1^2+1104*alpha__1-1133, index = real[2]), RootOf(1216*_Z^4+(264*alpha__2^8+408*alpha__2^7-1580*alpha__2^6-6832*alpha__2^5+3508*alpha__2^4+9944*alpha__2^3+9948*alpha__2^2-10752*alpha__2+5204)*_Z^3+(891*alpha__2^8+1652*alpha__2^7-4748*alpha__2^6-24076*alpha__2^5+5354*alpha__2^4+35356*alpha__2^3+29668*alpha__2^2-196*alpha__2+3971)*_Z^2+(506*alpha__2^8+980*alpha__2^7-2264*alpha__2^6-12420*alpha__2^5+3676*alpha__2^4+11596*alpha__2^3+33800*alpha__2^2-7772*alpha__2+1210)*_Z-473*alpha__2^8-720*alpha__2^7+2560*alpha__2^6+10960*alpha__2^5-8034*alpha__2^4-13840*alpha__2^3-9304*alpha__2^2+1104*alpha__2-1133, index = real[2]), RootOf(1216*_Z^4+(264*alpha__3^8+408*alpha__3^7-1580*alpha__3^6-6832*alpha__3^5+3508*alpha__3^4+9944*alpha__3^3+9948*alpha__3^2-10752*alpha__3+5204)*_Z^3+(891*alpha__3^8+1652*alpha__3^7-4748*alpha__3^6-24076*alpha__3^5+5354*alpha__3^4+35356*alpha__3^3+29668*alpha__3^2-196*alpha__3+3971)*_Z^2+(506*alpha__3^8+980*alpha__3^7-2264*alpha__3^6-12420*alpha__3^5+3676*alpha__3^4+11596*alpha__3^3+33800*alpha__3^2-7772*alpha__3+1210)*_Z-473*alpha__3^8-720*alpha__3^7+2560*alpha__3^6+10960*alpha__3^5-8034*alpha__3^4-13840*alpha__3^3-9304*alpha__3^2+1104*alpha__3-1133, index = real[2])}

 

{-0.7765721e-11, -0.40e-16, -0.2e-17}

(3)

evalf(bSol)

{.3768363389, .5974151794, 4.441922439}

(4)

Again, these may not be independent - this time Algfield fails

PA := evala(Algfield(bSol)); PA[4]

false

(5)

Could this be because of the real[2] indexing? Find the simple index corresponding to it. In each case it is index=1

evalf([seq(subs(real[2] = i, bSol[1]), i = 1 .. 4)]); evalf([seq(subs(real[2] = i, bSol[2]), i = 1 .. 4)]); evalf([seq(subs(real[2] = i, bSol[3]), i = 1 .. 4)])

[HFloat(0.3768363389232699), -HFloat(1.1781382393831348)+HFloat(1.6342731836523705)*I, HFloat(-0.8215019563741058), -HFloat(1.1781382393831348)-HFloat(1.6342731836523705)*I]

 

[HFloat(4.44192244093814), HFloat(0.7521012471815752)+HFloat(0.18644586398558147)*I, HFloat(-2.670902846485502), HFloat(0.7521012471815752)-HFloat(0.18644586398558147)*I]

 

[HFloat(0.5974151764626043), HFloat(8.500619958633818)+HFloat(3.0727707478333337)*I, HFloat(-3.6203669440591857), HFloat(8.500619958633818)-HFloat(3.0727707478333337)*I]

(6)

bSol2 := subs(real[2] = 1, bSol)

{RootOf(1216*_Z^4+(264*alpha__1^8+408*alpha__1^7-1580*alpha__1^6-6832*alpha__1^5+3508*alpha__1^4+9944*alpha__1^3+9948*alpha__1^2-10752*alpha__1+5204)*_Z^3+(891*alpha__1^8+1652*alpha__1^7-4748*alpha__1^6-24076*alpha__1^5+5354*alpha__1^4+35356*alpha__1^3+29668*alpha__1^2-196*alpha__1+3971)*_Z^2+(506*alpha__1^8+980*alpha__1^7-2264*alpha__1^6-12420*alpha__1^5+3676*alpha__1^4+11596*alpha__1^3+33800*alpha__1^2-7772*alpha__1+1210)*_Z-473*alpha__1^8-720*alpha__1^7+2560*alpha__1^6+10960*alpha__1^5-8034*alpha__1^4-13840*alpha__1^3-9304*alpha__1^2+1104*alpha__1-1133, index = 1), RootOf(1216*_Z^4+(264*alpha__2^8+408*alpha__2^7-1580*alpha__2^6-6832*alpha__2^5+3508*alpha__2^4+9944*alpha__2^3+9948*alpha__2^2-10752*alpha__2+5204)*_Z^3+(891*alpha__2^8+1652*alpha__2^7-4748*alpha__2^6-24076*alpha__2^5+5354*alpha__2^4+35356*alpha__2^3+29668*alpha__2^2-196*alpha__2+3971)*_Z^2+(506*alpha__2^8+980*alpha__2^7-2264*alpha__2^6-12420*alpha__2^5+3676*alpha__2^4+11596*alpha__2^3+33800*alpha__2^2-7772*alpha__2+1210)*_Z-473*alpha__2^8-720*alpha__2^7+2560*alpha__2^6+10960*alpha__2^5-8034*alpha__2^4-13840*alpha__2^3-9304*alpha__2^2+1104*alpha__2-1133, index = 1), RootOf(1216*_Z^4+(264*alpha__3^8+408*alpha__3^7-1580*alpha__3^6-6832*alpha__3^5+3508*alpha__3^4+9944*alpha__3^3+9948*alpha__3^2-10752*alpha__3+5204)*_Z^3+(891*alpha__3^8+1652*alpha__3^7-4748*alpha__3^6-24076*alpha__3^5+5354*alpha__3^4+35356*alpha__3^3+29668*alpha__3^2-196*alpha__3+3971)*_Z^2+(506*alpha__3^8+980*alpha__3^7-2264*alpha__3^6-12420*alpha__3^5+3676*alpha__3^4+11596*alpha__3^3+33800*alpha__3^2-7772*alpha__3+1210)*_Z-473*alpha__3^8-720*alpha__3^7+2560*alpha__3^6+10960*alpha__3^5-8034*alpha__3^4-13840*alpha__3^3-9304*alpha__3^2+1104*alpha__3-1133, index = 1)}

(7)

This time Algfield succeeds, and we only need the single algebraic number alpha__1

PA := evala(Algfield(bSol2)); PA[4]; K := PA[3]

true

 

{RootOf(11*_Z^9+17*_Z^8-64*_Z^7-280*_Z^6+142*_Z^5+370*_Z^4+376*_Z^3-96*_Z^2+47*_Z-11, index = 1)}

(8)

`~`[`@`(evala, Minpoly)](bSol2, _X)

{_X^9-(47/11)*_X^8+(96/11)*_X^7-(376/11)*_X^6-(370/11)*_X^5-(142/11)*_X^4+(280/11)*_X^3+(64/11)*_X^2-(17/11)*_X-1}

(9)

`~`[PolynomialTools[MinimalPolynomial]](bSol2, _X)

{11*_X^9-47*_X^8+96*_X^7-376*_X^6-370*_X^5-142*_X^4+280*_X^3+64*_X^2-17*_X-11}

(10)

NULL

Download minpoly.mw

It is a sum over the different roots. Perhaps there is an easier way, but here is one way to do it.

restart

S := proc (t) options operator, arrow; tanh(ln(1+t^2)) end proc

proc (t) options operator, arrow; tanh(ln(1+t^2)) end proc

q := diff(S(t), `$`(t, n))

pochhammer(1-n, n)+Sum(-(1/4)*_alpha^3*pochhammer(-n, n)*(t-_alpha)^(-1-n), _alpha = RootOf(_Z^4+2*_Z^2+2))

`assuming`([simplify(allvalues(q))], [n::posint])

(1/4)*GAMMA(n+1)*(-1)^n*((1+I)*(-1-I)^(1/2)*(t-(-1-I)^(1/2))^(-1-n)+(1-I)*(-1+I)^(1/2)*(t-(-1+I)^(1/2))^(-1-n)+(-1-I)*(-1-I)^(1/2)*(t+(-1-I)^(1/2))^(-1-n)+(-1+I)*(t+(-1+I)^(1/2))^(-1-n)*(-1+I)^(1/2))

``

Download simplify.mw

Your code contains the line

S2 := t -> piecewise(xn(t) <= xa and 0 < zn(t), -S1(t)*exp(mu1*alpha2(t)), 0)

which suggests that the drop-off might be when xn(t) or zn(t) becomes zero. Plotting xn(t) or equivalently, Xn, shows that Xn goes through zero at that point. So then fsolve(Xn, 0.2..0.9) returns the drop-off point of 0.5006805969.

min_problem.mw

Use Pi and not pi.

restart; with(LinearAlgebra)

w := c[i]*(1-cos(2*Pi*x/a))*(1-cos(2*Pi*y/b))

c[i]*(1-cos(2*Pi*x/a))*(1-cos(2*Pi*y/b))

al_eq1 := (1/2)*(int(int(D__11*(diff(w, x, x))^2+2*D__12*(diff(w, x, x))*(diff(w, y, y))+4*D__66*(diff(w, x, y))^2+D__22*(diff(w, y, y))^2-2*q__0*w, x = 0 .. a), y = 0 .. b))

c[i]*(6*D__11*Pi^4*b^5*c[i]+4*D__12*Pi^4*a^2*b^3*c[i]+6*D__22*Pi^4*a^4*b*c[i]+8*D__66*Pi^4*a^2*b^3*c[i]-a^4*b^5*q__0)/(a^3*b^4)

NULL

Download Pi.mw

A bit less efficient, perhaps.

[Edit: my interpretation of the OPs question was that different columns were to be scaled differently, but if rows were intended, then premultiplication by the diagonal matrix works.]

A:=Matrix(3,3,[1,6,9,7,4,3,2,8,9]);
B:=<a,b,c>;
A^~(-1).LinearAlgebra:-DiagonalMatrix(B);

Matrix(3, 3, {(1, 1) = 1, (1, 2) = 6, (1, 3) = 9, (2, 1) = 7, (2, 2) = 4, (2, 3) = 3, (3, 1) = 2, (3, 2) = 8, (3, 3) = 9})

Vector(3, {(1) = a, (2) = b, (3) = c})

Matrix(%id = 36893490033002534060)

``

Download div.mw

If I understand correctly what you want, a (Tabulated) DataFrame is one way to display this.

stiffnessmatrixcal.mw

First 30 31 32 33 34 35 36 Last Page 32 of 81