Advanced Search

dharr

Dr. David Harrington

8270 Reputation

22 Badges

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

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
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

not scanned...

dharr 8270
February 17 2019
0 3

I don't see a problem here - my interpretation is that seq returns the default value without accessing ("scanning") the entry in the rtable, just as V[6] returns 0.

exhaustive test...

dharr 8270
February 01 2019
2 4

I think the following logic works for mod 2 as well are regular matrices. We want to test all nonsingular matrices without too much duplication. Below shows that we need only test all perturbation matrices (n! of them) plus an extra 1 (n^2-n possibilities) and another extra 1 (n^2-n-1 possibilities). This is possible in about an hour on my laptop. If produces many duplicates and I show only one answer. It would be possible to eliminate duplicates with a bit more effort.

Of the 8 columns, none have only zeros so six of then have one 1 and two have two 1s. For the six columns with one 1, there will be two row positions, say i and j that have no 1s for any of the six columns. Now consider the 7th column, with two 1s.

If neither of these two 1s are in row i or j then this column is a sum of two of the six columns and so the matrix would be singular. So the 7th column has one or two 1's in row i or j. If it has only one, say in i, then for the 8th column there must be a 1 in row j, otherwise row j would be zero and the matrix would be singular. In this case, every column has a 1 in a different row (a permutation matrix) and there are two extras.

If the 1s in column 7 are in row i and j, then consider column 8. If its two 1s avoided row i and j, it would be a sum of the six columns. So it must have one 1 in i or j. Then column 7 has a 1 in the other or i or j and again we have a perturbation matrix with two extra 1s.

The following tests this.

> 

restart;
> 

with(LinearAlgebra[Modular]):with(Iterator):
> 

n:=8; #s = n+2
f:=x^8+x^7+x^5+x+1;#f:=x^8+x^4+x^3+x^2+1;
perms:=Permute(n):
Number(perms)*(n^2-n)*(n^2-n-1); # matrices to check

8

x^8+x^7+x^5+x+1

124185600

> 

st:=time():
p:=0:          # count matrices checked
nsucc:=0:      # count matrices with correct char. poly.
succ:=table(): # collect succesful ones
for perm in perms do
  A:=Create(2,n,n,integer[]):
  for i to n  do
    A[perm[i],i]:=1:
  end do:
  for k to n^2 do
    if A(k)=1 then next else A(k):=1 end if;
    for m to n^2 do
      if A(m)=1 then next else A(m):=1 end if;
      p:=p+1;
      #if p mod 10000000 =0 then print(p,(time()-st)/60*min) end if;
      if CharacteristicPolynomial(2,A,x)=f then nsucc:=nsucc+1; succ[nsucc]:=copy(A) end if;
      A(m):=0;
    end do;
    A(k):=0;
  end do;
end do:
print(p,`matrices tested.`,nsucc,`successes.`,`cpu time:`,(time()-st)/60*min);

124185600, `matrices tested.`, 80640, `successes.`, `cpu time:`, 67.80286667*min

> 

B:=succ[1];add(B);CharacteristicPolynomial(2,B,x);

_rtable[18446744628225754526]

10

x^8+x^7+x^5+x+1

> 

 


 

Download Modular5.mw

Edit: improved with modifications suggested by @vv :

Modular6.mw

parametric plots...

dharr 8270
January 29 2019
0 1

A bit contrived but this does it:


 

> 

A:=[1,2,3,4,5,6];
B:=[10,20,30,40,50,60];

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

[10, 20, 30, 40, 50, 60]

> 

plots:-display(seq(plot([A[i],y,y=0..B[i]],color=blue),i=1..numelems(A)),view=[0..6,DEFAULT],thickness=3);

> 

 


 

Download plots.mw

maybe not faster...

dharr 8270
January 29 2019
1 3

[Edit: see improved answer below]

Still takes a long time for 8x8 so I didn't wait. Below is just a different approach, FYI.
 

> 

restart;
> 

with(LinearAlgebra[Modular]):with(Iterator):
> 

n:=4;s:=3;
f:=x^4+x;
perms:=Permute([1$s,0$(n^2-s)]):
Number(perms);

4

3

x^4+x

560

> 

st:=time():i:=0:
for perm in perms do
 i:=i+1;
 lst:=convert(perm[],list);
 A:=Matrix(n,n,lst);
 cp:=CharacteristicPolynomial(2,A,x);
 if cp=f then print(time()-st,A,cp) end if;
end do:
i;
time()-st;

.281, Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (3, 1) = 0, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 1, (4, 4) = 0}), x^4+x

.328, Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 1, (4, 1) = 0, (4, 2) = 1, (4, 3) = 0, (4, 4) = 0}), x^4+x

.375, Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 1, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 1, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 1, (4, 4) = 0}), x^4+x

.406, Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 1, (2, 1) = 1, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 1, (4, 3) = 0, (4, 4) = 0}), x^4+x

.484, Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 1, (4, 1) = 1, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0}), x^4+x

.547, Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (2, 1) = 1, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0}), x^4+x

.578, Matrix(4, 4, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 1, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0}), x^4+x

.640, Matrix(%id = 18446744237432444678), x^4+x

560

.703

> 

 


 

Download Modular2.mw

Edit: since you wanted nonsingular I should have chosen x^4+x+1 (degree 4 without zero root) and s=4 (>=4 or there will be a zero row/column)

Edit: OP has changed algorithm so above is no longer comparable.

fsolve needs help...

dharr 8270
January 21 2019
0 0

I think you will have to give some ranges for the variables to make any progress here. These can probably be quite wide, e.g., 0..infinity if you know a parameter is positive and don't wan't to be too restrictive.

standard errors...

dharr 8270
January 21 2019
0 0

This is the first order approximation that @Mac Dude alludes to, with the details to get the standard errors out:
 

Download Parameters.mw

partly...

dharr 8270
January 16 2019
1 0

This goes partly there - need to collect the a's and convert to the equivalen binomial

> 

convert((a+b)^n,FormalPowerSeries,b) assuming n::posint;
convert(%,binomial);

Sum(a^n*(-1/a)^k*pochhammer(-n, k)*b^k/factorial(k), k = 0 .. infinity)

Sum(a^n*(-1/a)^k*binomial(-n+k-1, -n-1)*b^k, k = 0 .. infinity)

> 

 

> 

 


 

Download binomial.mw

Convert integrals to differential equati...

dharr 8270
January 06 2019
0 14

The problem is slow because at every step in z that the DE solver needs, it needs to evaluate many numerical integrals from the initial point to the current z position. I would suggest that you formulate the problem without the integrals, replacing them with DEs and then add these DEs to the existing two, and solve the larger system. For example, M1=Int(integrand with tt, tt=-500..z) can be diff(M1(z),z) = integrand with z with initial condition M1(-500)=0. Then the solver can evaluate M1(z) as well as the other quantities as it steps forward in z.

missing multiplication...

dharr 8270
December 28 2018
1 1

There is a missing multiplication. If this is inserted, then the indefinite integration is found, and Maple needs to know more (needs some assumptions made) to complete the definite integration.


 

> 

restart;
> 

integrand:=sin(x)/(a*b+a^2*sin(x)^2-d^2*cos(x)^2)*(c+cos(x));

sin(x)*(c+cos(x))/(a*b+a^2*sin(x)^2-d^2*cos(x)^2)

> 

indef:=int(integrand, x);

(1/2)*ln(a^2*cos(x)^2+d^2*cos(x)^2-a^2-a*b)/(a^2+d^2)-c*arctanh((a^2+d^2)*cos(x)/(a*(a+b)*(a^2+d^2))^(1/2))/(a*(a+b)*(a^2+d^2))^(1/2)

> 

int(integrand, x=0..X);

Warning, unable to determine if 2*Pi*_Z2+arccos((a*(a+b)*(a^2+d^2))^(1/2)/(a^2+d^2)) is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if 2*Pi*_Z3-arccos((a*(a+b)*(a^2+d^2))^(1/2)/(a^2+d^2)) is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if -Pi+2*Pi*_Z4+arccos((a*(a+b)*(a^2+d^2))^(1/2)/(a^2+d^2)) is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if 2*Pi*_Z2+Pi-arccos((a*(a+b)*(a^2+d^2))^(1/2)/(a^2+d^2)) is between 0 and X; try to use assumptions or use the AllSolutions option

int(sin(x)*(c+cos(x))/(a*b+a^2*sin(x)^2-d^2*cos(x)^2), x = 0 .. X)

> 

g:=simplify(eval(indef,x=X)-eval(indef,x=0));

(-2*c*(a^2+d^2)*arctanh((a^2+d^2)*cos(X)/(a*(a+b)*(a^2+d^2))^(1/2))+2*c*(a^2+d^2)*arctanh((a^2+d^2)/(a*(a+b)*(a^2+d^2))^(1/2))+(a*(a+b)*(a^2+d^2))^(1/2)*(ln((a^2+d^2)*cos(X)^2-a*(a+b))-ln(-a*b+d^2)))/((a*(a+b)*(a^2+d^2))^(1/2)*(2*a^2+2*d^2))

> 

simplify(g,symbolic);

-(1/2)*((a^2+d^2)^(1/2)*(ln(-a*b+d^2)-ln((a^2+d^2)*cos(X)^2-a*(a+b)))*a^(1/2)*(a+b)^(1/2)+2*(a^2+d^2)*(arctanh((a^2+d^2)^(1/2)*cos(X)/(a^(1/2)*(a+b)^(1/2)))-arctanh((a^2+d^2)^(1/2)/(a^(1/2)*(a+b)^(1/2))))*c)/(a^(1/2)*(a^2+d^2)^(3/2)*(a+b)^(1/2))

> 

 


 

Download integration.mw

2D math vs 2D input...

dharr 8270
December 21 2018
1 0

@chty_syq  You have several different "1+2" in this document, but only some are active. On the toolbar  "2D input" indicates that it looks nice but isn't active and "2D math" means that it is also active. "Maple input" means it is active but you are entering it in "1D" or command input style. You also want math rather than text on the line above.

If you open a new worksheet (rather than document) then enter 2-D math at the prompt, and press enter, the answer will appear and a new prompt for another math input, which is probably the easiest way to use Maple for math unless you want very sophisticated fomatting.

3j calculator...

dharr 8270
December 09 2018
1 0

The Maple applications centre has a worksheet to download that calculates 3j symbols:

https://www.maplesoft.com/applications/view.aspx?SID=4896

convert to cos...

dharr 8270
December 05 2018
1 0

If you want to get the cos form, then

int(tan(2*x),x);
convert(%,cos);
simplify(%,symbolic);

gives -(1/2)*ln(cos(2*x))

fsolve...

dharr 8270
December 05 2018
2 0

For example, to find up to 10 roots in the range -2..5, use

fsolve(eq,beta=-2..5,maxsols=10);

differentiate procedure...

dharr 8270
November 29 2018
2 2

You can just differentiate the procedure that produces Tf(x,y):

Tfnum:=rhs(op(3,pp));  # a procedure of y only

DTfnum:=D(Tfnum);

A fascinating capability of Maple.

DiffTf.mw

 

 

ODE output...

dharr 8270
November 22 2018
1 0

numSol, the output of dsolve can be used to get the numerical values out - try numSol(0.1). So extracting the values can be done something like this (depends on how exactly you want it formatted.

for tt from 0.1 to 0.9 by 0.1 do
  map(rhs,numSol(tt)[1..2])[];
end do;

First 65 66 67 68 69 70 71 Last Page 67 of 82