Carl Love

Carl Love

26488 Reputation

25 Badges

12 years, 299 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

The ssystem command returns a list of two items. The first item is an error code (0 meaning no error); the second item is a string containing the program's output. This string is the only part that you want.

The commands Import and ImportGraph have different syntax, and you are mixing them up. Import will take raw string data; ImportGraph requires a file or filename. Try this:

L2:= Import(
    ssystem("D:/nauty27r3/geng -c -b 6 -g")[2], 
    source= direct, format= "Graph6", output= list
);

Note that the command does not begin with "!".

However, you can verify that C6 is embeddable in S5 with a single command:

GT:= GroupTheory:
GT:-MinPermRepDegree(GT:-CyclicGroup(6));

                               5

 

Here's a procedure for it:
 

restart:

PolyToList:= (p::algebraic, V::list({name, function}))->
local T;
    sort(
        `[]`~([coeffs](collect(p, V, 'distributed'), V, T), map(degree~, [T], V)),
        #Sort by total degree with ties broken by left-to-right (wrt V) degrees:
        key= [add, op] @~ curry(op, 2)
    )
:

#Example usage:
V:= [x,y,sin(z)]:

P:= randpoly(V, terms= 16);

62*x^2-82*x*y+80*y^3-44*x^2*sin(z)^2+71*x*y^3-17*x*y*sin(z)^2-75*y^3*sin(z)-10*y*sin(z)^3-7*x^5-40*x^4*y+42*x^4*sin(z)-50*x^3*y*sin(z)+23*x^2*y*sin(z)^2+75*x^2*sin(z)^3-92*x*y^3*sin(z)+6*y*sin(z)^4

interface(rtablesize= 16):

PolyToList(P, V);

[[-82, [1, 1, 0]], [62, [2, 0, 0]], [80, [0, 3, 0]], [-10, [0, 1, 3]], [-75, [0, 3, 1]], [-17, [1, 1, 2]], [71, [1, 3, 0]], [-44, [2, 0, 2]], [6, [0, 1, 4]], [-92, [1, 3, 1]], [75, [2, 0, 3]], [23, [2, 1, 2]], [-50, [3, 1, 1]], [42, [4, 0, 1]], [-40, [4, 1, 0]], [-7, [5, 0, 0]]]

print~((Matrix@PolyToList)~(P, combinat:-powerset(V)[2..])):

Matrix(6, 2, {(1, 1) = 6*y*sin(z)^4-10*y*sin(z)^3-75*y^3*sin(z)+80*y^3, (1, 2) = [0], (2, 1) = -92*y^3*sin(z)-17*sin(z)^2*y+71*y^3-82*y, (2, 2) = [1], (3, 1) = 75*sin(z)^3+23*sin(z)^2*y-44*sin(z)^2+62, (3, 2) = [2], (4, 1) = -50*sin(z)*y, (4, 2) = [3], (5, 1) = 42*sin(z)-40*y, (5, 2) = [4], (6, 1) = -7, (6, 2) = [5]})

Matrix(3, 2, {(1, 1) = 75*x^2*sin(z)^3+42*x^4*sin(z)-7*x^5-44*x^2*sin(z)^2+62*x^2, (1, 2) = [0], (2, 1) = 6*sin(z)^4+23*x^2*sin(z)^2-50*sin(z)*x^3-40*x^4-10*sin(z)^3-17*sin(z)^2*x-82*x, (2, 2) = [1], (3, 1) = -92*sin(z)*x-75*sin(z)+71*x+80, (3, 2) = [3]})

Matrix(5, 2, {(1, 1) = -7*x^5-40*x^4*y+71*x*y^3+80*y^3+62*x^2-82*x*y, (1, 2) = [0], (2, 1) = 42*x^4-50*x^3*y-92*x*y^3-75*y^3, (2, 2) = [1], (3, 1) = 23*x^2*y-44*x^2-17*x*y, (3, 2) = [2], (4, 1) = 75*x^2-10*y, (4, 2) = [3], (5, 1) = 6*y, (5, 2) = [4]})

Matrix(10, 2, {(1, 1) = 6*sin(z)^4-10*sin(z)^3, (1, 2) = [0, 1], (2, 1) = -17*sin(z)^2-82, (2, 2) = [1, 1], (3, 1) = 75*sin(z)^3-44*sin(z)^2+62, (3, 2) = [2, 0], (4, 1) = -75*sin(z)+80, (4, 2) = [0, 3], (5, 1) = 23*sin(z)^2, (5, 2) = [2, 1], (6, 1) = -92*sin(z)+71, (6, 2) = [1, 3], (7, 1) = -50*sin(z), (7, 2) = [3, 1], (8, 1) = 42*sin(z), (8, 2) = [4, 0], (9, 1) = -40, (9, 2) = [4, 1], (10, 1) = -7, (10, 2) = [5, 0]})

Matrix(14, 2, {(1, 1) = 80*y^3, (1, 2) = [0, 0], (2, 1) = -75*y^3, (2, 2) = [0, 1], (3, 1) = 71*y^3-82*y, (3, 2) = [1, 0], (4, 1) = -92*y^3, (4, 2) = [1, 1], (5, 1) = 62, (5, 2) = [2, 0], (6, 1) = -10*y, (6, 2) = [0, 3], (7, 1) = -17*y, (7, 2) = [1, 2], (8, 1) = 6*y, (8, 2) = [0, 4], (9, 1) = 23*y-44, (9, 2) = [2, 2], (10, 1) = -50*y, (10, 2) = [3, 1], (11, 1) = -40*y, (11, 2) = [4, 0], (12, 1) = 75, (12, 2) = [2, 3], (13, 1) = 42, (13, 2) = [4, 1], (14, 1) = -7, (14, 2) = [5, 0]})

Matrix(11, 2, {(1, 1) = -7*x^5+62*x^2, (1, 2) = [0, 0], (2, 1) = 42*x^4, (2, 2) = [0, 1], (3, 1) = -40*x^4-82*x, (3, 2) = [1, 0], (4, 1) = -44*x^2, (4, 2) = [0, 2], (5, 1) = -50*x^3, (5, 2) = [1, 1], (6, 1) = 75*x^2, (6, 2) = [0, 3], (7, 1) = 23*x^2-17*x, (7, 2) = [1, 2], (8, 1) = 71*x+80, (8, 2) = [3, 0], (9, 1) = -10, (9, 2) = [1, 3], (10, 1) = -92*x-75, (10, 2) = [3, 1], (11, 1) = 6, (11, 2) = [1, 4]})

Matrix(%id = 18446746120599651622)

 

Download PolyToList.mw

I hope this helps you:
 

NULL

restart:

f:= x-> c*min(x, 1-x):

#
#Find fixed points of f:
f1x:= {solve}(x=f(x), x);

{0, c/(c+1)}

#Find 2-orbits of f, which you observed as .4 and .8:
f2x:= {solve}(x=f(f(x)), x);

{0, c/(c+1), c/(c^2+1), c^2/(c^2+1)}

f2x:= f2x minus f1x:  f1x:= [f1x[]]:  f2x:= [f2x[]]:

c:= 2:

evalf[1](f2x);

[.4, .8]

plot(
    [f@@~[0$2,1,2][], ([`$`]~)~([f1x,f2x], 2)[]], 0..1,
    style= [line$4, point$2], symbol= [solidcircle, soliddiamond],
    color= [white, red, blue, green, red, blue], symbolsize= 20,
    legend= [
        ``, #just for blank line at top of legend
        typeset~(f^~``~([0,1,2]), "\n")[],
        typeset~("\n", ["fixed point", "2-orbit"], " of  ", f, " \n")[]
    ],
    legendstyle= [location= right],
    axes= boxed, view= [(-.1..1.1)$2], scaling= constrained, size= [600$2],
    tickmarks= [(f@@0=typeset)~([f1x[], f2x[], 1])$2], gridlines
);

NULL

NULL


 

Download OrbitPlot.mw

Like this:

z:= [
    0, 0, 0, 0, 0, 0, .689376362, 1.378752724, 2.068129087, 2.757505449, 0, 1.02920355,
    2.0584071, 3.087610649, 4.116814199, 0, 1.216469264, 2.432938529, 3.649407793,
    4.865877057, 0, 1.325720912, 2.651441823, 3.977162735, 5.302883646
]:
plots:-listcontplot(
    [ListTools:-LengthSplit](z,5), axis[1,2]= [tickmarks= ([$1..5]=~[$0..4])]
);

 

The Maple and Mathematica solutions that you show are mathematically equivalent. Both will give nonreal values for any real values of the a's and b's (unless a[1,1]*b[0,2] = a[0,2]*b[1,1], in which case the discriminant is 0). The expression under the square root sign cannot be positive for any real values of the parameters because it factors as (-1)*(a[1,1]*b[0,2] - a[0,2]*b[1,1])^2.

The absence of from an expression is not a reliable indicator that the expression is real-valued. The presence of is not a reliable indicator that it is not real-valued.

For rational-coefficient polynomial equations (such as yours) with parameters, you can get solve to break down all the relevant real intervals of the parameters which produce real solutions by adding the options parametric and real:

solve(
    k^2*a[1,1]^2 + k^2*b[1,1]^2 + 4*k*a[0,2]*a[1,1] + 
    4*k*b[0,2]*b[1,1] + 4*a[0,2]^2 + 4*b[0, 2]^2 = 0, 
    k, parametric, real
);

Note that the solution below covers all degenerate cases, i.e., when the equation is not really quadratic because the coefficient of k^2 is 0 and even when the coefficient of k is also 0. The parametric option always considers those cases.
   

         

I agree that a list or sequence of vectors would be the preferred format for programmatic usage. For a displayable table, I prefer using a matrix. (I especially like the way that matricies display with their entries centered in both dimensions.) The code below does both: My procedure DiffTab returns the "raw" list of vectors, the same as Tom Leslie's doDiff. The procedures PadV and DisplayTab turn those vectors into a neatly displayable matrix.

PadV:= proc(V::Vector, n::posint)
local k:= numelems(V), OP:= j-> n-k-1+2*j;
    Vector(2*n-1, applyop~(OP, 1, op(2,V)), fill= ``)
end proc
:
DisplayTab:= (DT::list(Vector), X::Vector)-> `<|>`(PadV~([X, DT[]], numelems(X))[])
:
DiffTab:= proc(Y::Vector)
local k, r:= Y, R:= table([1=r]);
    for k from 2 to numelems(Y) do r:= r[2..] - r[..-2];  R[k]:= r od;
    convert(R, list)
end proc
:
X:= <[$1..5]>:
Y:= 1/~X:
DT:= DiffTab(Y);
<<` X` | ` Y` | [``$4][]>, DisplayTab(DT,X)>;

             

You need to change xy to x*y.

Note the diffferent output from these two:

restart:
solve(y(1)^2 = 2., y(1));

             
1.414213562, -1.414213562

{ fsolve(y(1)^2 = 2., y(1)) };
             { }

The fsolve returns no output (which I emphasized by putting the command in { }), not even an error message. Technically, the fsolve output is called NULL.

An expression like y(1) is called in Maple an unevaluated function call (or sometimes just a "function"). So, you can use a function as a variable for solve but not for fsolve. So, your call T__v_x(1000) is producing NULL output in a position where some other command is expecting a number.

 

Do you mean something like this?

#x and y coordinates of grid:
GX:= [seq](0..3, 1/3):  GY:= [seq](1..3, 1/3)
:
#x and y coordinates of highlighted points:
PX:= [seq](1/4..11/4, 1/4):  PY:= [seq](5/4..11/4, 1/4)
:
plots:-display(
    #border:
    plot(
        [[GX[1],GY[1]], [GX[1],GY[-1]], [GX[-1],GY[-1]], [GX[-1],GY[1]], [GX[1],GY[1]]],    
        color= COLOR(RGB, .3$3), thickness= 2
    ),
    #inner:
    plot(
        [
            seq([[GX[1],y],[GX[-1],y]], y= GY[2..-2]), 
            seq([[x,GY[1]],[x,GY[-1]]], x= GX[2..-2])
        ],
        color= COLOR(RGB, .5$3), thickness= 0
    ),
    #points:
    plot(
        [seq](seq([x,y], x= PX), y= PY),
        style= point, symbol= solidcircle, symbolsize= 8,
        color= red
    ),
    axes= none, scaling= constrained
);

Your general solution isn't general enough. If f is *any* differentiable function such that f(0)=1, then x0 = y0 = f is a solution to your equation. 

A basic fact about systems of equations is that the number of things that can be solved for is at most the number of equations. 

If is the list of expressions, then all that you need to do is

R:= Grid:-Map(int, L, x= 0..1);

For numeric integration, do 

R:= Grid:-Map(int, L, x= 0..1, numeric);

The command will automatically apply some rudimentary load balancing.

It's possible to get the exact solutions for all 6 functions by direct use (i.e., without using dsolve) of the Method of Undetermined Coefficients. If you try to do this with dsolve, it will get the first 5, but the 6th took longer than I cared to wait. I killed it after more than an hour. The method below gets it in seconds.
 

Direct application of Method of Undetermined Coefficients

restart
:

#right sides of the 6 odes:
RHS:= [
    0, -v__0^3, -3*v__0^2*v__1, -3*v__0*v__1^2 - 3*v__0^2*v__2,
    -6*v__0*v__1*v__2 - 3*v__0^2*v__3 - v__1^3,
    -6*v__0*v__1*v__3 - 3*v__1^2*v__2 - 3*v__0^2*v__4 - 3*v__0*v__2^2
]:
print~(RHS)
:

0

-v__0^3

-3*v__0^2*v__1

-3*v__0^2*v__2-3*v__0*v__1^2

-3*v__0^2*v__3-6*v__0*v__1*v__2-v__1^3

-3*v__0^2*v__4-6*v__0*v__1*v__3-3*v__0*v__2^2-3*v__1^2*v__2

ICs:= [[1,0], [0,0]$5] #initial conditions of the 6 odes
:

for k to nops(RHS) do
    S:= C1*sin(t) + C2*cos(t); #general homogenous solution
    #Perform method of undetermined coefficients termwise:
    rs:= combine(expand(RHS[k]));
    for f in indets(rs, specfunc({sin,cos})) do #for each term on right side
        C:= coeff(rs, f);
        d:= degree(C,t);
        o:= op(f);
        P:= add(a[i]*t^i, i= 0..d)*sin(o) + add(b[i]*t^i, i= 0..d)*cos(o);
        if o=t then P*= t fi; #extra power of t   
        S+= eval(P, solve(identity(diff(P,t$2)+P=C*f, t)))
    od;
    #Apply initial conditions to determine C1 and C2:
    v__||(k-1):= (combine@expand@eval)(S, solve(eval([S, diff(S,t)], t=0)=~ICs[k]));
    print(v[k-1] = v__||(k-1))
od:

v[0] = cos(t)

v[1] = -(1/32)*cos(t)-(3/8)*sin(t)*t+(1/32)*cos(3*t)

v[2] = (23/1024)*cos(t)+(3/32)*sin(t)*t-(3/128)*cos(3*t)+(1/1024)*cos(5*t)-(9/128)*cos(t)*t^2-(9/256)*t*sin(3*t)

v[3] = -(547/32768)*cos(t)+(9/1024)*sin(t)*t^3-(207/4096)*sin(t)*t+(135/4096)*cos(t)*t^2+(279/8192)*t*sin(3*t)-(81/4096)*t^2*cos(3*t)+(297/16384)*cos(3*t)-(3/2048)*cos(5*t)+(1/32768)*cos(7*t)-(15/8192)*t*sin(5*t)

v[4] = (1/1048576)*cos(9*t)-(99/16384)*sin(t)*t^3-(9/131072)*cos(7*t)+(1539/65536)*t^2*cos(3*t)+(825/262144)*t*sin(5*t)-(1359/65536)*cos(t)*t^2-(3915/131072)*t*sin(3*t)+(883/524288)*cos(5*t)-(15121/1048576)*cos(3*t)+(8997/262144)*sin(t)*t+(6713/524288)*cos(t)+(27/32768)*cos(t)*t^4-(225/131072)*t^2*cos(5*t)+(243/32768)*t^3*sin(3*t)-(21/262144)*t*sin(7*t)

v[5] = -(3/1048576)*cos(9*t)+(4635/1048576)*sin(t)*t^3-(81/1310720)*sin(t)*t^5+(1757/16777216)*cos(7*t)-(96795/4194304)*t^2*cos(3*t)-(16575/4194304)*t*sin(5*t)+(15777/1048576)*cos(t)*t^2+(108243/4194304)*t*sin(3*t)-(921/524288)*cos(5*t)+(394701/33554432)*cos(3*t)-(108357/4194304)*sin(t)*t-(42397/4194304)*cos(t)-(441/4194304)*t^2*cos(7*t)-(27/8388608)*t*sin(9*t)+(2187/1048576)*t^4*cos(3*t)+(1125/1048576)*t^3*sin(5*t)+(1/33554432)*cos(11*t)-(783/1048576)*cos(t)*t^4+(6975/2097152)*t^2*cos(5*t)-(10935/1048576)*t^3*sin(3*t)+(1659/8388608)*t*sin(7*t)

 

Download UndeterminedCoeffs.mw

The plot command is plotting over a real interval. It can't be expected to find on its own every turn of a function, especially considering that there could be infinitely many of them (such as for the Weierstrass function[*1]). If there are values of n that you want to make sure are included, you can use the sample option. In this case, we might as well make the sample all the integers in the interval:

f:= n-> ceil(sqrt(4*floor(n))) - floor(sqrt(2*floor(n))) - 1:

Edit: My intention was to use the OP's original function, but I mistakenly copy and pasted Tom Leslie's modification of the original. So, make that

f:= n-> ceil(sqrt(4*n)) - floor(sqrt(2*n)) - 1:
plot(f, 10...100, sample= [$10..100]);

[*1] Weierstrass function: f:= x-> Sum(cos(3^n*x)/2^n, n= 0..infinity). It is continuous everywhere and differentiable nowhere.
 

You should use for derivatives in initial and boundary conditions. I don't know why that is, but it was the only change that I needed to make to your code. Sometimes the diff form works, but the form always works and is a lot less to type.

sys:= {
    diff(phi(eta), eta$2) + 5.261282735*f(eta)*diff(phi(eta), eta) - 
    2.630641368*phi(eta) = 0,
 
    1.059704409*diff(theta(eta), eta$2) + 6.176017503*f(eta)*diff(theta(eta), eta) 
    + 21.03607964*diff(f(eta), eta$2) + 0.5*phi(eta) = 0,
 
    diff(f(eta), eta$4) - 1.052256547*diff(f(eta), eta)*diff(f(eta), eta$2) + 
    1.052256547*f(eta)*diff(f(eta), eta$3) + 5.165076420*diff(theta(eta), eta) + 
    5.261282735*diff(phi(eta), eta) = 0,

    D(phi)(0) = 1 + 0.5*(D@@2)(f)(0),
    D(phi)(1) = 0.5*(D@@2)(f)(1), 
    f(0) = -0.5, f(1) = 0.5, 
    phi(0) = 1, phi(1) = 0, 
    theta(0) = 1, theta(1) = 0
}:
dsol:= dsolve(sys, numeric):
plots:-odeplot(
    dsol, `[]`~(eta, [f, phi, theta](eta)), eta= 0..1, 
    legend= [f, phi, theta]
);

First 35 36 37 38 39 40 41 Last Page 37 of 382