Advanced Search

vv

13827 Reputation

20 Badges

9 years, 317 days
Contact vv

MaplePrimes Activity


These are answers submitted by vv

restart...

vv 13827
November 25 2017
0 2

Put restart in a separate execution group.

You may also switch to Assume using

with(MathematicalFunctions,Assume);
Assume(s>0);

 

Domains...

vv 13827
November 25 2017
0 2

For working with polynomials over GF, the simplest method is to use the Domains package .

This package does not seem to be very popular, and even some links in the documentation are broken., e.g. the link Domains/example should be
?examples,Domains.

 

If you are interested only in the field Fp i.e. GF(p,1) then Domains is not necessary because GF(p,1) is Zp and the operator mod is enough.

 

For example, to compute the gcd of two polynomials in Z7[x]:

 

> 

a := x^5+5*x^4+6*x+2;
b := x^5+3*x^3+5*x^2+2*x+5;

x^5+5*x^4+6*x+2

  

x^5+3*x^3+5*x^2+2*x+5

 (1)
> 

Gcd(a,b) mod 7;

x^2+1

 (2)

Note that the gcd of these polynomials in Z[x] is 1

> 

gcd(a,b);

1

 (3)

Now let's compute the same gcd using Domains:

> 

restart;
> 

with(Domains);

---------------------- Domains version 1.0 ---------------------
Initially defined domains are Z and Q the integers and rationals
Abbreviations, e.g. DUP for DenseUnivariatePolynomial, also made

  

[AbelianGroup, AbelianMonoid, AlgebraicExtension, Array, CommutativeRing, CramersRule, DenseExponentVector, DenseUnivariatePolynomial, DifferentialField, DistributedMultivariatePolynomial, EuclideanAlgorithm, EuclideanDomain, ExpandedNormalFormQuotientField, ExponentVector, FactoredNormalFormQuotientField, Field, Finite, FiniteField, FrobeniusNormalForm, GaloisField, Gaussian, GaussianInteger, GcdDomain, Generator, GrobnerBasis, Group, Integers, IntegralDomain, LazyUnivariatePowerSeries, List2Vector, MacaulayExponentVector, Maple, MapleDistributedMultivariatePolynomial, MapleExponentVector, MapleGaussianArithmetic, MapleMultivariatePolynomial, MapleUnivariatePolynomial, Matrix, MatrixDomain, Modp1Polynomial, Modp2Polynomial, Monoid, New, NormalForm, ODESolve, OrderedAbelianGroup, OrderedAbelianMonoid, OrderedDomain, OrderedField, OrderedSet, OrderedUnivariatePolynomial, PartiallyOrderedSet, PolynomialHomomorphism, Poset, PowerRemainder, PrimeExponentVector, PrincipalIdeal, Q, QuotientField, RandMatrix, RationalFunction, Rationals, RelativelyPrimeHeuristic, RepeatedSquaring, Ring, SemiGroup, Set, SingleExtensionGaloisField, Singular, Solve, SparseDistributedMultivariatePolynomial, SparseMultivariatePolynomial, SquareFreeNorm, SquareMatrix, TableDistributedMultivariatePolynomial, TranscendentalFunctions, UniqueFactorizationDomain, UnivariatePolynomial, UnivariatePowerSeries, Vector, Z, Zmod, addCategory, addProperties, addProperty, defOperation, defOperations, evaldomains, hasCategories, hasCategory, hasOperation, hasOperations, hasProperties, hasProperty, isDomain, isImplemented, newCategory, newDomain, notImplemented, show]

 (4)
> 

C := DenseUnivariatePolynomial(GF(7,1), x):
> 

a := x^5+5*x^4+6*x+2;

x^5+5*x^4+6*x+2

 (5)
> 

b := x^5+3*x^3+5*x^2+2*x+5;

x^5+3*x^3+5*x^2+2*x+5

 (6)
> 

a1:=C[Input](a);

`domains/DenseUnivariatePolynomial/badge0`(2, 6, 0, 0, 5, 1)

 (7)
> 

b1:=C[Input](b);

`domains/DenseUnivariatePolynomial/badge0`(5, 2, 5, 3, 0, 1)

 (8)
> 

C[Gcd](a1,b1);

`domains/DenseUnivariatePolynomial/badge0`(1, 0, 1)

 (9)

Note trhat the conversion to a1 and b1 is necessary:

> 

C[Gcd](a,b);

`domains/DenseUnivariatePolynomial/badge0`(1)

 (10)
> 

lprint(a);

x^5+5*x^4+6*x+2

  
> 

lprint(a1);

`domains/DenseUnivariatePolynomial/badge0`(2, 6, 0, 0, 5, 1)

  

 

A last remark. If you want to implement your own gcd, then note that the elements of GF(p,k) are themselves polynomials (actually in the quotient ring Zp[alpha]/(r), where r is an irreducible polynomial in the variable alpha, of degree k).  So, those  zppolys are elements of GF, not the polynomials in x].

 

By change of variables...

vv 13827
November 23 2017
0 7
> 

rel:=simplify(tan((x-t)/(alpha*(t^2+1)))-t);

-tan((-x+t)/(alpha*(t^2+1)))-t

 (1)
> 

alpha:=1/2;

1/2

 (2)
> 

X:=solve(rel,x);

(1/2)*t^2*arctan(t)+t+(1/2)*arctan(t)

 (3)
> 

t1:=fsolve(X=1);

.6177384029

 (4)
> 

evalf(Int(1/(1+t^2)*diff(X,t),t=0..t1));

.8903158470

 (5)
> 

 

Edit: without a change (i.e. directly):

restart;
rel:=simplify(tan((x-t)/(alpha*(t^2+1)))-t):
alpha:=1/2:
T:=solve(rel,t):
Int(1/(1+T^2),x=0..1):
evalf(%);

       .8903158471

First two eigenvalues...

vv 13827
November 22 2017
0 2 
R:=i -> RootOf(_Z^9-28*_Z^8+328*_Z^7-2088*_Z^6+(-2*cos(y)-2*cos(x)+7852)*_Z^5+(28*cos(y)+28*cos(x)-17768)*_Z^4+(-152*cos(y)-152*cos(x)+23632)*_Z^3+(408*cos(y)+408*cos(x)-17232)*_Z^2+(-12*cos(x)*cos(y)-540*cos(y)-540*cos(x)+5844)*_Z+32*cos(x)*cos(y)+256*cos(y)+256*cos(x)-544, index=i):
plot3d([seq(R(i),i=1..2)], x=0..2*Pi,y=0..2*Pi);

unapply...

vv 13827
November 22 2017
1 0

For the second situation you must use

f__c:=unapply(piecewise(cnd, 1, not cnd, 0),  x__0, y__0);

That's due to the evaluation rules for procedures. But in essence, x__0, y__0 in cnd are global variables, distinct from those in proc(x__0,y__0); they are just "homonymes".

polylog(1/2, -1) + 1...

vv 13827
November 21 2017
3 7

The sum is computed correctly to polylog(1/2, -1) + 1
Instead of stopping here, Maple tries to use
    polylog(a, z) + polylog(a, -z) = 2^(1-a)*polylog(a, z^2)
for z=1 ==> division by 0.

 

global t...

vv 13827
November 20 2017
1 2

In the DynamicSystems package Sinc() depends on the global variable t. So, you must use t for its value. For t=0, Sinc is not defined and generates an error.

with(DynamicSystems):
Sinc();
      sin(Pi*t)/(Pi*t)

eval(Sinc(),t=1);
    0

eval(Sinc(),t=0);
Error, numeric exception: division by zero

limit(Sinc(), t=0);
    1

eval(Sinc(), t=0.0001);
    .9999999834

 

 

 

 

add...

vv 13827
November 19 2017
0 4

The approximation of the sum of a series is an almost impossible task.
So, I can understand that Maple is counting on the user's ability to use add (and assume the responsibility!).

Otherwise, we obtain:

evalf( Sum(piecewise(n<29,1/n^2,1/n), n=1..infinity) );
    unevaluted

evalf( Sum(piecewise(n<30,1/n^2,1/n), n=1..infinity) );
    1.644934067

 

 

 

VectorCalculus...

vv 13827
November 18 2017
1 4

Actually you have computed  J(F)' . Grad(G)
where J demotes the Jacobian and ' is the transpose.
You may use the VectorCalculus package which works for n variables and any coordinate system:

 

> 

restart;
> 

with(VectorCalculus):
> 

J:=Jacobian([f(x,y,z,t),g(x,y,z,t),h(x,y,z,t),k(x,y,z,t)], [x,y,z,t]);

_rtable[18446744074328080502]

 (1)
> 

dG:=Matrix(Gradient(G(x,y,z,t),[x,y,z,t]));

_rtable[18446744074369828734]

 (2)
> 

J^+ . dG;

_rtable[18446744074369819342]

 (3)
> 

 


Download J.mw

solpart...

vv 13827
November 18 2017
0 5

Here is a procedure (not tested enough).
The system must be polynomial with rational coefficients.

solpart:=proc(sys::list(relation))
local i,r1,r2,rez:=NULL, sd:=SolveTools:-SemiAlgebraic(sys);
#print(sd);
if sd=[] then return [] else sd:=sd[1] fi;
for i to nops(sd) do
   r1:=eval(sd[i],[rez]); 
   if type(r1,`=`) then rez:=rez,r1; next; fi;
   if i=nops(sd) then
     if     type(rhs(r1),realcons) then rez:=rez,lhs(r1)=rhs(r1)-1
     elif   type(lhs(r1),realcons) then rez:=rez,rhs(r1)=lhs(r1)+1
     else error "Strange1"; fi;
   return rez;
   fi; 
   
   r2:=eval(sd[i+1],[rez]);   #print([r1,r2]);
   if    type(rhs(r1),name) and (rhs(r1)=lhs(r2)) then i:=i+1; rez:=rez, rhs(r1)=eval((lhs(r1)+rhs(r2))/2,[rez]) 
   elif  type(lhs(r1),name) and (lhs(r1)=rhs(r2)) then i:=i+1; rez:=rez, lhs(r1)=eval((rhs(r1)+lhs(r2))/2,[rez])
   elif  type(rhs(r1),realcons) then rez:=rez,lhs(r1)=rhs(r1)-1
   elif  type(lhs(r1),realcons) then rez:=rez,rhs(r1)=lhs(r1)+1
   else error "Strange2"; fi;
od;
rez
end:

solpart([x*y-1<=0, x>0, y > 1]);
        y = 2, x = 1/4

solpart([x*y-1<0, x>0, y >= 1]);
        y = 1, x = 1/2

 

Why matrix?...

vv 13827
November 16 2017
1 0 
f:=proc(n::posint,i)
   local k, s:=add(i[k],k=1..n), p:=mul(i[k],k=1..n);
   add(p/i[k]*(s-i[k]), k=1..n)
end:

f(4,i);
      

InertForm...

vv 13827
November 16 2017
1 1

In Maple it is difficult to control the order of the terms in a sum (the sort command is mainly for polynomials).  
Maple does this for efficiency reasons; after all it was designed for computations rather than typesetting (where LaTeX is anyway the preferred tool by mathematicians).

In recent versions it offers however many typesetting options (see ?Typesetting).

 

Another method would be to use InertForm, but I am not sure whether it exists/works in Maple 13.

 

> 

restart;
> 

p:=A/(s+1)+B/(s-2)+C/(s+3);

A/(s+1)+B/(s-2)+C/(s+3)

 (1)
> 

with(InertForm):
> 

ip:=MakeInert(p):
> 

ip1:=subs([A=4/15,B=8/15,C=-3/10],ip);

`%+`(`%/`(4/15, `%+`(s, 1)), `%/`(8/15, `%+`(s, -2)), `%/`(-3/10, `%+`(s, 3)))

 (2)
> 

Display(ip=ip1);

`%+`(`%/`(A, `%+`(s, 1)), `%/`(B, `%+`(s, -2)), `%/`(C, `%+`(s, 3))) = `%+`(`%/`(4/15, `%+`(s, 1)), `%/`(8/15, `%+`(s, -2)), `%/`(-3/10, `%+`(s, 3)))

 (3)

The correct value...

vv 13827
November 13 2017
0 0 
Int(1/(sqrt(t^4+1)+t^2), t = 1 .. infinity):
% = value(%);
evalf(%); #check

        0.4799537050 = 0.4799537050

2^18...

vv 13827
November 11 2017
1 1

Such a system will have generally 2^18 = 262144 solutions.
E.g. the simple one
{seq(x[k]^2=1,k=1..18)};

has exactly 2^18 solutions (all are real). For a more realistic system each solution will be huge.
So, you have no chance. You can at most use fsolve to find some of them.

 

 

 

solve...

vv 13827
November 08 2017
1 1

solve(convert(eq-K,list));

First 80 81 82 83 84 85 86 Last Page 82 of 120