Advanced Search

dharr

Dr. David Harrington

8492 Reputation

22 Badges

21 years, 34 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

randomly multiply by matrices...

dharr 8492
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 8492
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 8492
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 8492
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 8492
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 8492
July 01 2020
0 0

You need a more accurate calculation - see the attached
 

Download Digits.mw

convert to D...

dharr 8492
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 8492
June 29 2020
1 0

Another way:

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

 

partly using the geometry package...

dharr 8492
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 8492
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 8492
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?

too many ICs...

dharr 8492
June 24 2020
0 1

Your ode_y differential equation is second order in time, so can have only two initial conditions, but you supplied 12.

a guess...

dharr 8492
June 23 2020
1 0

I'm using Maple 2017, but I don't think that matters. I could probably check my guess below by looking at the code, but that's a lot of work.
 

> 

restart;
> 

infolevel[dsolve]:=5;

5

 (1)
> 

ode:=diff(y(x),x)=x*ln(y(x));
> 

dsolve(ode,y(x));

diff(y(x), x) = x*ln(y(x))

  

Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature
trying 1st order linear
trying Bernoulli
trying separable
<- separable successful

  

y(x) = exp(RootOf(x^2+2*Ei(1, -_Z)+2*_C1))

 (2)
> 

dsolve([ode,y(1)=1],y(x)); #no user information printed

y(x) = 1

 (3)
> 

dsolve([ode,y(1)=3],y(x)); #but there is here

 -> Computing symmetries using: way = 3

 -> Computing symmetries using: way = 2
 -> Computing symmetries using: way = patterns
   -> The symmetry found is [1/x 0]
Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature
trying 1st order linear
trying Bernoulli
trying separable
<- separable successful

  

y(x) = 1/exp(RootOf(x^2+2*Ei(1, _Z)-2*Ei(1, -ln(3))-1))

 (4)

So I'm guessing that it tries the special case y(x)=y(1) and never gets to a deeper analysis.


 

Download ode.mw

currentdir after saving...

dharr 8492
June 21 2020
1 1

There can be some confusion about currentdir() as follows. If I open Maple from a shortcut/start menu etc and in a blank worksheet type currentdir() then it tells me some location where my Maple installation is; probably not a place I want to save a file. If I then save the worksheet somewhere useful, say on my desktop, close Maple and then open the worksheet I've just saved by double-clicking on it, now currentdir() shows my desktop location and files will be saved there.

convert menu...

dharr 8492
June 20 2020
0 0

I tend to type 1D math into 2D math input regions, and so something like Int(x^2,x=0..infinity) (where I typed right-arrow after the 2) produces the right output, but the Int and infinity don't change to their nice forms. If I want to change them as I go, then I use ^ (or @acer's escape is better), but if I want to change the whole thing to look better, I use the right-click menu 2-D math/convert to/2D math input, and then everything will be converted to the nice form. (Or I start with a 1D math execution group, do it all in 1D math, then convert it). So convert will work after you paste 1D math into a 2D region.

I don't know of any way to convert all instances in a worksheet.

First 57 58 59 60 61 62 63 Last Page 59 of 83