Advanced Search

dharr

Dr. David Harrington

8210 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
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

with(plots)...

dharr 8210
July 11 2020
2 1

You need with(plots): since implicitplot, pointplot and display are all in the plots package. Using with(plots,implicitplot) loads only implicitplot and not the others.

0 or 1...

dharr 8210
July 11 2020
2 0

about(_B1) returns
 Originally _B1, renamed _B1~:
  is assumed to be: OrProp(0,1)

so it means either zero or 1

D...

dharr 8210
July 10 2020
2 0


 

> 

restart;

> 

with(plots): with(VectorCalculus):
rf:=t-> < cos(t), sin(t), t >;

proc (t) options operator, arrow; VectorCalculus:-`<,>`(cos(t), sin(t), t) end proc

> 

rCurve := spacecurve( rf(t), t = -1 .. 1 );

> 

drf:=D(rf):
> 

drCurve := spacecurve( drf(t), t = -1 .. 1 );

> 

 

 

Download D.mw

addcoords only works one way...

dharr 8210
July 09 2020
3 0

to use addcoordinates you need to know x, y and z in terms of the new coordinate system, it doesn't do inverse ones. pdechangecoords is deprecated, but dchange can do what you want, though the form of the inverse transformation is a mess:
 

> 

restart;
> 

with(PDEtools):
> 

toroidaltrans:={x=a*sinh(v)*cos(w)/(cosh(v)-cos(u)),y=a*sinh(v)*sin(w)/(cosh(v)-cos(u)),z=a*sin(u)/(cosh(v)-cos(u))};

{x = a*sinh(v)*cos(w)/(cosh(v)-cos(u)), y = a*sinh(v)*sin(w)/(cosh(v)-cos(u)), z = a*sin(u)/(cosh(v)-cos(u))}

Solving tor the toroidal oordinates in terms of the cartesian coordinates gives a mess. Adding assumptions and using simplify with options might help here.

> 

invtoroidaltrans:=solve(toroidaltrans,{u,v,w});

{u = arctan(-2*a*z/((2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2+2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)*RootOf((-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2+a^2+x^2+y^2+z^2)*_Z^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)), (a^2-x^2-y^2-z^2)/((2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2+2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)*RootOf((-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2+a^2+x^2+y^2+z^2)*_Z^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2))), v = ln(RootOf((-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2+a^2+x^2+y^2+z^2)*_Z^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)), w = arctan(y*RootOf((x^2+y^2)*_Z^2-1), RootOf((x^2+y^2)*_Z^2-1)*x)}

> 

fn:=diff(f(x,y,z),z);

diff(f(x, y, z), z)

To toroidal

> 

fntoroidal:=dchange(toroidaltrans,fn,[u,v,w],params={a});

(-(diff(f(u, v, w), u))+(diff(f(u, v, w), u))*cosh(v)*cos(u)-(diff(f(u, v, w), v))*sinh(v)*sin(u))/a

And back again

> 

fnback:=dchange(invtoroidaltrans,fntoroidal,[x,y,z],itr=toroidaltrans,params={a},simplify);

diff(f(x, y, z), z)

> 

 


 

Download toroidal.mw

randomly multiply by matrices...

dharr 8210
July 07 2020
2 0

Perhaps something like this...

> 

restart;
> 

R[0]:=Matrix(3,3,[1,0,0,0,cos(delta),sin(delta),0,-sin(delta),cos(delta)]);
R[1]:=Matrix(3,3,[cos(delta),0,-sin(delta),0,1,0,sin(delta),0,cos(delta)]);

R[0] := Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = cos(delta), (2, 3) = sin(delta), (3, 1) = 0, (3, 2) = -sin(delta), (3, 3) = cos(delta)})

R[1] := Matrix(3, 3, {(1, 1) = cos(delta), (1, 2) = 0, (1, 3) = -sin(delta), (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = sin(delta), (3, 2) = 0, (3, 3) = cos(delta)})

Random number generators to select which matrix and which delta value

> 

r01:=rand(0..1):r2Pi:=rand(0..2.*Pi):

Successively multiply by a randomly selected Matrix with a randomly selected delta

> 

N:=1000;
S:=table():S[1]:=Vector([0,1,0]); #initial vector
for i from 2 to N do S[i]:=eval(R[r01()],delta=r2Pi()).S[i-1] end do:
S:=convert(S,list):

N := 1000

S[1] := Vector(3, {(1) = 0, (2) = 1, (3) = 0})

> 

plots:-pointplot3d(S,symbolsize=10,color=black);

> 

 


Edit: altered to make initial vector a unit vector.

Download Stokes.mw

partial answer...

dharr 8210
July 07 2020
0 1

I came to the same conclusion as @vv, that there is not a unique solution. Here is a solution to the simpler problem in which P6 is absent and length L25 is known.

> 

restart;

Mapleprimes problem proposed by @Fancypants.

Set up some viable points just to figure out some viable lengths. Make a drawing of the problem.
Problem: the lengths of the lines are given, as are P0 and P1. Find the locations of the other points.

> 

with(geometry):
pxy:={p3x=-0.5,p3y=3,p5x=3,p5y=4,p6x=3.5,p6y=3}:
point(P0,0,0):point(P1,1,0):point(P3,eval([p3x,p3y],pxy)):point(P5,eval([p5x,p5y],pxy)):point(P6,eval([p6x,p6y],pxy)):
OnSegment(P2,P0,P3,1.2):OnSegment(P4,P3,P5,0.7):
segment(P0P1,P0,P1):segment(P0P3,P0,P3):segment(P3P5,P3,P5):segment(P5P6,P5,P6):segment(P1P4,P1,P4):segment(P2P6,P2,P6):
plots:-display(draw([P0,P1,P2,P3,P4,P5,P6],printtext=true,axes=none),
               draw([P0P1,P0P3,P3P5,P5P6,P1P4,P2P6],axes=none));
lengths=[L01=distance(P0,P1),L02=distance(P0,P2),L23=distance(P2,P3),L34=distance(P3,P4),L45=distance(P4,P5),L56=distance(P5,P6),L26=distance(P2,P6),L14=distance(P1,P4),L25=distance(P2,P5)];

lengths = [L01 = 1, L02 = 1.658935235, L23 = 1.382446030, L34 = 1.498846154, L45 = 2.141208791, L56 = 1.118033989, L26 = 4.011605067, L14 = 3.412271768, L25 = 4.037018784]

Now solve it only knowing these lengths.

> 

restart;
> 

with(geometry):
> 

assume(L01>0,L02>0,L23>0,L34>0,L45>0,L14>0,L56>0,L26>0,L25>0,p3y>0,p5y>0,p6y>0);

Given lengths

> 

lengths:={L01 = 1, L02 = 1.658935235, L23 = 1.382446030, L34 = 1.498846154, L45 = 2.141208791, L56 = 1.118033989, L26 = 4.011605067, L14 = 3.412271768, L25 = 4.037018784}:
#assign(lengths);

P0 and P1 are known, set P0 at the origin and P1 at (L01,0)

> 

point(P0,0,0);point(P1,L01,0);

P0

P1

The points we need to find

> 

point(P3,p3x,p3y);point(P5,p5x,p5y);point(P6,p6x,p6y);

P3

P5

P6

Now define some segments we need to draw for the two quadrilaterals (the segments are just for drawing, not for calculation).

> 

segment(P0P1,P0,P1);segment(P0P3,P0,P3);segment(P3P5,P3,P5);segment(P5P6,P5,P6);

P0P1

P0P3

P3P5

P5P6

Point P2 is on line P0P3 at ratio k = L02/L23 and similarly for point P4

> 

OnSegment(P2,P0,P3,L02/L23);
OnSegment(P4,P3,P5,L34/L45);

P2

P4

Now we can define the other two segments

> 

segment(P1P4,P1,P4);segment(P2P6,P2,P6);

P1P4

P2P6

Setup the equations to solve for the locations of all the points. P0 and P1 are known immediately in terms of L01, and P2 and P4 are given in terms of the others, so  the 6 unknowns are the coordinates of P3, P5 and P6.

Specifying all the lengths gives only 5 independent equations - need another

> 

eq02:=L02^2=simplify(distance(P0,P2)^2);
eq23:=L23^2=simplify(distance(P2,P3)^2); #equiv to eq02
eq34:=L34^2=simplify(distance(P3,P4)^2);
eq45:=L45^2=simplify(distance(P4,P5)^2); #equiv to eq34
eq14:=L14^2=simplify(distance(P1,P4)^2);
eq56:=L56^2=simplify(distance(P5,P6)^2);
eq26:=L26^2=simplify(distance(P2,P6)^2);
eq25:=L25^2=simplify(distance(P2,P5)^2); #not supposed to know this length

L02^2 = L02^2*(p3x^2+p3y^2)/(L23+L02)^2

L23^2 = L23^2*(p3x^2+p3y^2)/(L23+L02)^2

L34^2 = (p3y^2-2*p3y*p5y+p5y^2+(p3x-p5x)^2)*L34^2/(L45+L34)^2

L45^2 = L45^2*(p3y^2-2*p3y*p5y+p5y^2+(p3x-p5x)^2)/(L45+L34)^2

L14^2 = ((L01^2-2*L01*p5x+p5x^2+p5y^2)*L34^2+2*L45*(L01^2+(-p3x-p5x)*L01+p3x*p5x+p3y*p5y)*L34+L45^2*(L01^2-2*L01*p3x+p3x^2+p3y^2))/(L45+L34)^2

L56^2 = p5y^2-2*p5y*p6y+p6y^2+(p5x-p6x)^2

L26^2 = (L02^2*(p3x^2-2*p3x*p6x+p3y^2-2*p3y*p6y+p6x^2+p6y^2)-2*(-p6y^2+p3y*p6y+p6x*(p3x-p6x))*L23*L02+L23^2*(p6x^2+p6y^2))/(L23+L02)^2

L25^2 = (L02^2*(p3x^2-2*p3x*p5x+p3y^2-2*p3y*p5y+p5x^2+p5y^2)-2*(-p5y^2+p3y*p5y+p5x*(p3x-p5x))*L23*L02+L23^2*(p5x^2+p5y^2))/(L23+L02)^2

First solve the slightly simpler problem, where we know the length L25 but don't have point P6 in the problem. Now we have 4 equations and 4 unkowns

> 

eqs:={eq02,eq34,eq14,eq25}:nops(eqs);
indets(eqs,name);

4

{L01, L02, L14, L23, L25, L34, L45, p3x, p3y, p5x, p5y}

Solve takes forever, so try fsolve, which returns the correct solution

> 

fsolve(eval(eqs,lengths),{p3x=-2..0,p3y=0..10,p5x=0..10,p5y=0..10});

{p3x = -.4999999995, p3y = 3.000000000, p5x = 2.999999999, p5y = 4.000000002}

Is there a unique solution with the information given, i.e. is the system rigid?

Adding the point P6 and replacing P2P5 by two other sides of the triangle P2P5P6 seems to allow for some flexibility?

 

Download geometry2.mw

longest distance with the path...

dharr 8210
July 05 2020
1 7

Slower than my earlier version but does find the path. The labelling method finds the edges in the path easily, but putting them in sequence to make a trail is not so simple.

> 

restart;

P.L. Giscard and P. Rochet, Graphs and Combinatorics, 34 (2018) 1197-1202, doi: 10.1007/s00373-018-1966-9 ; Eq 2 with W=zA where A is labelled and only subsets with a and b are used.

> 

LongestPath:=proc(a,b,G)
description "Find longest path between a and b in graph G";
local sumz,In,vertices,subvertices,invertices,subg,subgraph,H,A,N,j,subsets,zA,W,omega,plength,P,p,edges;
uses GraphTheory;
vertices:=Vertices(G);
N:=numelems(vertices);
W:=Matrix(N,N,symbol=omega);
subvertices:={vertices[]} minus {a,b};
subsets:=combinat:-subsets(subvertices):
sumz:=0:In:=LinearAlgebra:-IdentityMatrix(N):
while not subsets[finished] do
  invertices:=subsets[nextvalue]() union {a,b};
  subg:=InducedSubgraph(G,invertices);
  subgraph:=Graph(vertices,Edges(subg));
  H:=numelems(invertices);    
  A:=AdjacencyMatrix(subgraph)*~W;
  zA:=z*A;
  sumz:=sumz+zA^(H-1).(In-zA)^(N-H);
end do:
P:=collect(expand(sumz[a,b]),z); #Generating polynomial for number of paths of different lengths between a and b in graph G
plength:=degree(P,z);
p:=lcoeff(P,z);
# now find edges of longest path
if type(p,`+`) then p:=op(1,p) end if; #first one if more than 1
if not type(p,`*`) then
    edges:={{op(1..2,p)}}
  else
    edges:=map(x->{op(1..2,x)},{op(p)});
  end if;
plength,edges #return length of longest path and set of edges of a longest path
end proc:
> 

with(GraphTheory):
> 

G := Graph([1, 2, 3, 4, 5], {{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 1}, {1, 4}});

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(%id = 18446744782322561022), `GRAPHLN/table/1`, 0)

Longest path distance and edges in a longest path

> 

d,edges:=LongestPath(2,4,G);
HighlightEdges(G,edges);

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

> 

DrawGraph(G);

> 

G2:=RandomGraphs:-RandomGraph(10,16);

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], Array(%id = 18446744782322587998), `GRAPHLN/table/19`, 0)

> 

d,edges:=LongestPath(2,4,G2);
HighlightEdges(G2,edges);

8, {{1, 3}, {1, 6}, {2, 9}, {3, 5}, {4, 6}, {5, 8}, {8, 10}, {9, 10}}

> 

DrawGraph(G2);

> 

 

> 

 

 

Download PathPolyijlabelled.mw

longest distance but not the path...

dharr 8210
July 04 2020
3 1

Much slower than @vvs, and doesn't find the path, but does find all the distances.

> 

restart;

P.L. Giscard and P. Rochet, Graphs and Combinatorics, 34 (2018) 1197-1202, doi: 10.1007/s00373-018-1966-9 ; Eq 2 with W=zA and only subsets with i and j (a and b)

> 

PathPoly:=proc(a,b,G)
description "Generating polynomial for number of paths of different lengths between a and b in graph G";
local sumz,In,vertices,subvertices,invertices,subg,subgraph,H,A,N,j,subsets,zA;
uses GraphTheory;
vertices:=Vertices(G);
N:=numelems(vertices);
subvertices:={vertices[]} minus {a,b}; #all subsets with a and b
subsets:=combinat:-subsets(subvertices):
sumz:=0:In:=LinearAlgebra:-IdentityMatrix(N):
while not subsets[finished] do
  invertices:=subsets[nextvalue]() union {a,b};
  subg:=InducedSubgraph(G,invertices);
  subgraph:=Graph(vertices,Edges(subg));
  H:=numelems(invertices);    
  A:=AdjacencyMatrix(subgraph);
  zA:=z*A;
  sumz:=sumz+zA^(H-1).(In-zA)^(N-H);
end do:
collect(sumz[a,b],z);
end proc:
> 

with(GraphTheory):
> 

G := Graph([1, 2, 3, 4, 5], {{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 1}, {1, 4}});

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5], Array(%id = 18446744831140302710), `GRAPHLN/table/1`, 0)

> 

DrawGraph(G);

Two paths of length 2 and one of length 3. Longest path distance is the degree

> 

P:=PathPoly(2,4,G);
degree(P);

z^3+2*z^2

3

> 

G2:=RandomGraphs:-RandomGraph(10,16);

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], Array(%id = 18446744831179633718), `GRAPHLN/table/19`, 0)

> 

P2:=PathPoly(2,4,G2);
degree(P2);
totalpaths:=eval(P2,z=1);

2*z^8+4*z^7+6*z^6+4*z^5+z^4+z^3+2*z^2+z

8

21

> 

DrawGraph(G2);

> 

 

 

Download LongestPath4.mw

reading text files and include...

dharr 8210
July 01 2020
1 0

ModuleProcs.txt

TestModule.txt

TestModule.txt contains:

 

Test:=module()
  option package;
  export Lth;
  local Sqr;
$include "ModuleProcs.txt"
end module;

 

ModuleProcs.txt contains:

 

Sqr:=proc(A) A^2 end proc;
Lth:=proc(A,B,usersqrt) usersqrt(Sqr(A)+Sqr(B)) end proc;

 

> 

restart;
> 

read "TestModule.txt"

_m948865488512

Save the library

> 

savelib('Test',cat(currentdir(),"/Testlib.mla"));

Sqr is used by Lth without the user needing to know about it. The following is a user function to be used by Lth

> 

mysqrt:=proc(x) sqrt(x) end proc;

proc (x) sqrt(x) end proc

> 

Test:-Lth(3,5,mysqrt);

34^(1/2)

> 

 


 

Download module2.mw

Digits...

dharr 8210
July 01 2020
0 0

You need a more accurate calculation - see the attached
 

Download Digits.mw

convert to D...

dharr 8210
June 30 2020
2 0 
convert(ode,D);

gives

(D(y))(x) = x^2+x*y(x)

you can convert back by convert( ,diff)

This seems not quite what you expect. D(y)(x) means the differential of the function y evaluated at x. Perhaps you wanted something else.

Edit: D[](y)(x); is y(x). (differentiate zero times)

subsop...

dharr 8210
June 29 2020
1 0

Another way:

expr:=A*B*C;
subsop(1=1,expr);

 

partly using the geometry package...

dharr 8210
June 27 2020
1 3

The geometry package is good on various intersection points etc, but I used @Carl Love's chordal segment formula to finish off.

> 

restart;
> 

with(geometry):

Let's make the origin at A, with R=1 fixing the scale of the problem

> 

assume(L>0,L<2);
#L:=1.158728473; #uncomment for plotting purposes

Scale so A is at the origin and R=1. cA and cC are circles centered at A and C respectively.

> 

point(A,0,0);point(C,-1,0);
circle(cA,[A,1]); # center A, radius 1
circle(cC,[C,L]); # center C, radius L

A

C

cA

cC

Intersection points of the two circles are G and H, and we want the distance between them

> 

intersection('GandH',cC,cA,[G,H]);
line(lGH,[G,H]);
dGH:=distance(G,H);

[G, H]

lGH

(4-(-L^2+2)^2)^(1/2)

Maple isn't so good with areas, so for the area of the circular segment will use @Carl Love's routine

> 

AreaChord:= proc(L,R) local x;2*int(sqrt(R^2-x^2) - sqrt(R^2-(L/2)^2), x= 0..L/2) end:

The area of the two segments is required to be half the area of circle cA

> 

eq:=AreaChord(dGH,1)+AreaChord(dGH,L)=area(cA)/2;
LL:=fsolve(eq,L);

-(1/2)*L*(-L^2+4)^(1/2)+arcsin((1/2)*L*(-L^2+4)^(1/2))+L^2*arcsin((1/2)*(-L^2+4)^(1/2)) = (1/2)*Pi

1.158728473

To see the diagram rerun with the correct L value at the top of the worksheet

> 

plots:-display(
  draw([A,C,G,H],printtext=true),
  draw([cA,cC,lGH]),
axes=none);

> 

 


 

Download circles.mw

stretched circle...

dharr 8210
June 27 2020
0 0

Adding to @Carl Love's method, you can just plot a circle in a square, and then use the view option of plot to have different ranges for the vertical and horizontal axes.

helper procedure...

dharr 8210
June 24 2020
0 2

Not sure I fully understand. But if you just want Dist to access Sqr then

Dist := proc (X) Sqr(X) end proc;

and then

Test:-Dist(b)

will return b^2. Is this what you want?

First 55 56 57 58 59 60 61 Last Page 57 of 81