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

matrix and graph generation...

dharr 8270
April 05 2017
0 0

I would do the generation of the matrices and graphs like this (for the 2x2 case)

> 

n:=2;
bits:=n^2;
mlist:=table():
for i from 0 to 2^bits-1 do
   mlist[i]:=Matrix(n,n,Bits:-Split(i,'bits'=bits));
end do:
mlist:=convert(mlist,list);

n := 2

  

bits := 4

  

mlist := [Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 0, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 0, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 1, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 1, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 1, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 1, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 1, (2, 2) = 1})]

 (1)
> 

graphs:=map(GraphTheory:-Digraph,mlist):
> 

 


 

Download graphs.mw

indexed procedures...

dharr 8270
March 04 2017
0 1

You can set up indices with procedures, but you stll need parentheses to show it is a procedure.

> 

q := proc () local inds; if not type(procname, 'indexed') then return 'procname(args)' end if; inds := op(procname); inds end proc;

proc () local inds; if not type(procname, 'indexed') then return 'procname(args)' end if; inds := op(procname); inds end proc

 (1)
> 

q(y)

q(y)

 (2)
> 

q[4, 5]()

4, 5

 (3)
> 

q[4, 5]

q[4, 5]

 (4)
> 

``

 

Flawed logic...

dharr 8270
February 17 2017
0 20

If r=0 then you correctly change i. But then you test r agoinst N+1, which fails, and then the else clause sets i back to r, which is zero.

Loop weight...

dharr 8270
February 17 2017
0 1

The Digraph command doesn't seem to allow loops, but adding a weight in a diagonal entry of the Adjacency matrix seems to be accepted, and the weight is given by DrawGraph near the vertex.


 

> 

with(GraphTheory):
> 

A := Matrix([[2,2,0],[0,0,0],[1,3,0]]);

A := Matrix(3, 3, {(1, 1) = 2, (1, 2) = 2, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 1, (3, 2) = 3, (3, 3) = 0})

 (1)
> 

G := Graph(A);

GRAPHLN(directed, weighted, [1, 2, 3], Array(%id = 18446744074240237686), `GRAPHLN/table/20`, Matrix(%id = 18446744074240237566))

 (2)
> 

DrawGraph(G);

 
> 

Digraph({[1,1]});

Error, (in GraphTheory:-Graph) invalid edge/arc: [1, 1]

  
> 

 


 

Download Loopgraph.mw

Make some assumptions...

dharr 8270
February 11 2017
0 1

If you add "assuming M>0, g>0;" then you get undefined. If your further choose a limit from the left or right you get -infinity or +infinity.

Give the list of the new variables...

dharr 8270
February 11 2017
0 16

Maple needs to know the new varables are u and v, not W and Q. Try:

Change(V, {x = u-W, y = v-Q}, [u, v]);

Another way...

dharr 8270
December 06 2016
0 2

if Row(A,2).Row(A,2)=0 then...

missing multiplication...

dharr 8270
December 05 2016
0 2

Another 2-D math problem. There needs to be a multiplication between the c^2 and the parenthesis. Then it easily plots, and fsolve gives the two solutions, which are obviously x=+/-R.

y(x)2.mw

Not unique...

dharr 8270
November 30 2016
1 4

The companion matrix can be used - for monic polynomials. The determinant is not unique and may not be the one you want.
 

> 

with(LinearAlgebra):
> 

p:=w^2 + u*w + v; #monic polynomial

u*w+w^2+v

 (1)
> 

w*IdentityMatrix(degree(p,w))-CompanionMatrix(p,w);

Matrix([[w, v], [-1, w+u]])

 (2)
> 

 


 

Download Companion.mw

There is no exact solution, only a close...

dharr 8270
November 02 2016
2 2

There is no exact solution to the accuracy that Maple usually works at. You need to find a close "near enough" solution, which means you can't use fsolve. The easiest way is to minimize a sum of squares as in the worksheet below.
 

> 

eqs := {(2^2/a^2+2^2/c^2-2*(2*2)*cos(x)/(a*c))/sin(x)^2 = 1/2.2393^2, (2^2/a^2+4^2/c^2+4*(2*2)*cos(x)/(a*c))/sin(x)^2 = 1/1.5968^2, (1^2/a^2+sin(x)^2/b^2+2^2/c^2+(2*2)*cos(x)/(a*c))/sin(x)^2 = 1/2.7896^2, 1/sin(x)^2*(2^2*sin(x)^2/b^2) = 1/2.8650^2};

{(4/a^2+4/c^2-8*cos(x)/(a*c))/sin(x)^2 = .1994230893, (4/a^2+16/c^2+16*cos(x)/(a*c))/sin(x)^2 = .3921922000, (1/a^2+sin(x)^2/b^2+4/c^2+4*cos(x)/(a*c))/sin(x)^2 = .1285038476, 4/b^2 = .1218290191}

 (1)

Move the sin(x)^2 over to the right hand side

> 

eqs2:=eqs*~sin(x)^2;

{4*sin(x)^2/b^2 = .1218290191*sin(x)^2, 4/a^2+4/c^2-8*cos(x)/(a*c) = .1994230893*sin(x)^2, 4/a^2+16/c^2+16*cos(x)/(a*c) = .3921922000*sin(x)^2, 1/a^2+sin(x)^2/b^2+4/c^2+4*cos(x)/(a*c) = .1285038476*sin(x)^2}

 (2)

Sunm of squares of differences to minimize

> 

S:=add((lhs(eq)-rhs(eq))^2,eq in eqs2);

(4*sin(x)^2/b^2-.1218290191*sin(x)^2)^2+(4/a^2+4/c^2-8*cos(x)/(a*c)-.1994230893*sin(x)^2)^2+(4/a^2+16/c^2+16*cos(x)/(a*c)-.3921922000*sin(x)^2)^2+(1/a^2+sin(x)^2/b^2+4/c^2+4*cos(x)/(a*c)-.1285038476*sin(x)^2)^2

 (3)

Minimise subject to some constraints

> 

Optimization:-Minimize(S,{a>=0,b>=0,c>=0,x>=0,x<=Pi});

[0.162642134569434096e-11, [a = HFloat(8.902766366460824), b = HFloat(5.7300076115735505), c = HFloat(6.416618564610261), x = HFloat(1.8419029174006967)]]

 (4)

Gives a close solution, but not the one you wanted (though b is right). Probably there are several close solutions; it will be hard to decide - perhaps there are several combinations of a and c and x that work

> 

 


 

Download Crystal.mw

fsolve only returns one answer...

dharr 8270
November 01 2016
1 2

In general fsolve only returns one answer. If you want others, you can specify a range to look for roots in, e..g, x=187.5..188.5 or give an initial guess, e.g, x=188. See the help for fsolve

Loop rather than fsolve...

dharr 8270
October 25 2016
2 2

Using fsolve to find the number of terms at which the error is 0.3 won't work because the error for an integer number of terms won't exactly be 0.3, so there is no solution. You know that the error must decrease so it is efficient to just add an extra term on and then recalculate the error in a loop.

The integration method _d01akc works well for these sorts of integrals.

> 

restart;
> 

f:=x->x;
y:=(n,x)->exp(-x/2)*sin(n*Pi*x);
w:=(x)->exp(x);

proc (x) options operator, arrow; x end proc

  

proc (n, x) options operator, arrow; exp(-(1/2)*x)*sin(n*Pi*x) end proc

  

proc (x) options operator, arrow; exp(x) end proc

 (1)

Note the normalizing integral is always 1/2:

> 

d:=int(w(x)*y(n,x)^2,x=0..1) assuming n::posint;

1/2

 (2)

And we can get a general formula for the coefficients

> 

int(w(x)*f(x)*y(n,x),x=0..1)/d assuming n::posint;
c:=unapply(%,n);

-8*n*Pi*(4*Pi^2*exp(1/2)*(-1)^n*n^2+3*(-1)^(1+n)*exp(1/2)+4)/(16*Pi^4*n^4+8*Pi^2*n^2+1)

  

proc (n) options operator, arrow; -8*n*Pi*(4*Pi^2*exp(1/2)*(-1)^n*n^2+3*(-1)^(n+1)*exp(1/2)+4)/(16*Pi^4*n^4+8*Pi^2*n^2+1) end proc

 (3)
> 

Fourierf:=0:Lerror:=infinity:
for n while Lerror>0.3 do
  Fourierf:=Fourierf+c(n)*y(n,x);
  Lerror:=evalf(sqrt(Int((f(x)-Fourierf)^2*w(x),x=0..1,method='_d01akc'))); #edit: test after adding term
  print(n,Lerror);
od:

1, .5894736138

  

2, .4627928802

  

3, .3943241015

  

4, .3483437147

  

5, .3156572254

  

6, .2905059507

 (4)
> 

 

 

Download Eigenfn_error.mw

Works in Maple 2015...

dharr 8270
September 06 2016
0 0

In Maple 2015.2 evalf(Int(...)) quickly returns the answer you give.

expand within sin...

dharr 8270
July 28 2016
2 4


> 

q:=sin(omega*(T0-T)+Phi);

sin(omega*(T0-T)+Phi)

 (1)
> 

map(expand,q);algsubs(T0*omega+Phi=c,%);expand(%);subs(c=T0*omega+Phi,%);

sin(-T*omega+T0*omega+Phi)

  

-sin(T*omega-c)

  

-sin(omega*T)*cos(c)+cos(omega*T)*sin(c)

  

-sin(omega*T)*cos(T0*omega+Phi)+cos(omega*T)*sin(T0*omega+Phi)

 (2)
> 

 

 

 

 

Bits package...

dharr 8270
April 02 2016
1 2

Bits:-GetBits(19,-1..0,bits=8);

0, 0, 0, 1, 0, 0, 1, 1

Bits:-GetBits(-(19+1),-1..0,bits=8);

1, 1, 1, 0, 1, 1, 0, 0  #ones complement

Bits:-GetBits(-19,-1..0,bits=8);

1, 1, 1, 0, 1, 1, 0, 1  #twos complement

First 69 70 71 72 73 74 75 Last Page 71 of 82