Advanced Search

dharr

Dr. David Harrington

8507 Reputation

22 Badges

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

Equivalent Resistance...

dharr 8507
February 28 2022
2 0

I started off doing this with the correct number of node and loop equations, using CycleBasis for the loop equations. But by defining a voltage per node rather than voltage per edge, the loop equations can be skipped. There are more equations than unknowns, but solve can handle the redundancy.

Resistance:=proc(G::Graph,VerticesInOut::set)
  local dG,VertexIn,VertexOut,nV,nE,currents,voltages,
    currInOut,nodeeqns,deviceeqns,arcs,k,i,v,x,ans;
  uses GraphTheory;
  VertexIn,VertexOut:=VerticesInOut[];
  arcs:=map(convert,Edges(G),list);
  dG:=Graph(Vertices(G),arcs); #directed version required for signed incidence matrix
  # Define voltages per node. Current arrow direction for edge {j,k} is j->k (j<k)
  nV:=NumberOfVertices(G);
  nE:=NumberOfEdges(G);
  currents:=[seq(i[j],j=1..nE)];
  voltages:=[seq(v[j],j=1..nV)];
  # unit current in and out through VertexIn and VertexOut 
  currInOut:=Vector(nV,0);
  currInOut[VertexIn]:=1;
  currInOut[VertexOut]:=-1;
  # Node equations from incidence matrix - rows are vertices, columns are edges.
  # One of these equations is redundant, but solve can handle that
  nodeeqns:={LinearAlgebra:-GenerateEquations(IncidenceMatrix(dG,reverse),currents,currInOut)[]};
  # Device equations - Ohms law with R=1
  deviceeqns:=table();
  for k to nops(arcs) do:
    x:=arcs[k];
    deviceeqns[k]:=v[x[1]]-v[x[2]]=i[k];
  end do;
  deviceeqns:=convert(deviceeqns,set);
  # Solve the equations
  ans:=solve(nodeeqns union deviceeqns,{currents[],voltages[]});
  eval(v[VertexIn]-v[VertexOut],ans)
end proc:
> 

G:=GraphTheory:-SpecialGraphs:-SoccerBallGraph();

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], Array(%id = 18446744073861270334), `GRAPHLN/table/1`, 0)

> 

Resistance(G,{1,31});
Resistance(G,{1,2});

17/11

16273/25080

> 

 

 

 

resistanceZ.mw

assume...

dharr 8507
February 19 2022
0 0

Maybe assuming it is a constant is what you want.

> 

D(d*x);

D(d)*x+d*D(x)

> 

`assuming`([D(d*x)], [d::constant]);

d*D(x)

> 

assume(a, constant);
> 

D(a*x);

a*D(x)

> 

 

Download constant.mw

same engine...

dharr 8507
February 18 2022
2 0

You can use the same Maple engine for both worksheets, and then the Matrix stays in memory from the first to the second. Under Tools->Options in the general tab under mathematical engine choose Share one engine among all documents. Of course if you have variables with the same name in both worksheets they have to have the same purpose.

Maple-assisted...

dharr 8507
February 12 2022
0 0

Don't know much about this, and didn't try to make a procedure, just did it step by step with Maple.

NULL

See https://mathworld.wolfram.com/Sigma-Algebra.html. for definition of sigma algebra

> 

restart
> 

X := {1, 2, 3}

{1, 2, 3}

Suppose C is our developing collection, which starts with:

> 

C := {{1}, {2}}

{{1}, {2}}

Complements must be in C

> 

map(proc (i) options operator, arrow; `minus`(X, i) end proc, C)
C := `union`(C, %)

{{1, 3}, {2, 3}}

{{1}, {2}, {1, 3}, {2, 3}}

For any sequence in C its union must be in C - assume null set counts as a sequence

> 

seqs := combinat:-powerset(C)

{{}, {{1}}, {{2}}, {{1, 3}}, {{2, 3}}, {{1}, {2}}, {{1}, {1, 3}}, {{1}, {2, 3}}, {{2}, {1, 3}}, {{2}, {2, 3}}, {{1, 3}, {2, 3}}, {{1}, {2}, {1, 3}}, {{1}, {2}, {2, 3}}, {{1}, {1, 3}, {2, 3}}, {{2}, {1, 3}, {2, 3}}, {{1}, {2}, {1, 3}, {2, 3}}}

All unions of sequences must be in C

> 

sequnions := map(proc (x) options operator, arrow; `union`(x[]) end proc, seqs)
C := `union`(C, sequnions)

{{}, {1}, {2}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}

{{}, {1}, {2}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}

Maybe there are some complements of these we need to add? yes, {3} is new

> 

newcomps := map(proc (i) options operator, arrow; `minus`(X, i) end proc, C)
C := `union`(C, newcomps)

{{}, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}

{{}, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}

We have now all subsets so there can't be more

> 

nops(C) = nops(combinat:-powerset(X))

8 = 8

> 

NULL

NULL

 

Download sigma_algebra.mw

Maple formatting for Notepad++...

dharr 8507
February 12 2022
2 1

Here are two links. I haven't used either so can't comment further.

https://www.mapleprimes.com/posts/129552-maple-in-all-colours-of-rainbow-

https://github.com/yanglr/HighlightForMaple

dsolve...

dharr 8507
February 09 2022
3 0

I'm not exactly sure if I've understood. I ignored any dc part, since it didn't seem to be present in your plots (though I wasn't clear what U(t) was supposed to be). Assuming you want to solve analytically and not numerically, here is my take on it.
 

> 

restart

Series R-L-C

> 

params := {C = 1/401, L = 100, R = 20, V__0 = 2000}

{C = 1/401, L = 100, R = 20, V__0 = 2000}

Device equations

> 

eqL := V__L(t) = L*(diff(q(t), t, t));

V__L(t) = L*(diff(diff(q(t), t), t))

q(t) = C*V__C(t)

V__R(t) = (diff(q(t), t))*R

Kirchoff

> 

eq := V__0 = V__R(t)+V__C(t)+V__L(t);

V__0 = V__R(t)+V__C(t)+V__L(t)

Initial condtions

> 

ics := (D(q))(0) = V__0/L, q(0) = 0;

(D(q))(0) = V__0/L, q(0) = 0

> 

ans := dsolve({eq, eqC, eqL, eqR, ics}, {V__C(t), V__L(t), V__R(t), q(t)});

{V__C(t) = (-(1/2)*exp(-(1/2)*(R*C-(C*(C*R^2-4*L))^(1/2))*t/(C*L))*(R^2*C^2+C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L-2*(C^2*R^2-4*C*L)^(1/2))*V__0/(C*R^2-4*L)-(1/2)*exp(-(1/2)*(R*C+(C*(C*R^2-4*L))^(1/2))*t/(C*L))*V__0*(R^2*C^2-C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L+2*(C^2*R^2-4*C*L)^(1/2))/(C*R^2-4*L)+V__0*C)/C, V__L(t) = -(1/2)*((1/2)*R^2/L+(1/2)*(C*(C*R^2-4*L))^(1/2)*R/(C*L)-1/C)*V__0*(R^2*C^2-C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L+2*(C^2*R^2-4*C*L)^(1/2))*exp(-(1/2)*(R*C+(C*(C*R^2-4*L))^(1/2))*t/(C*L))/(C*R^2-4*L)-(1/2)*((1/2)*R^2/L-(1/2)*(C*(C*R^2-4*L))^(1/2)*R/(C*L)-1/C)*(R^2*C^2+C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L-2*(C^2*R^2-4*C*L)^(1/2))*V__0*exp(-(1/2)*(R*C-(C*(C*R^2-4*L))^(1/2))*t/(C*L))/(C*R^2-4*L), V__R(t) = -(1/2)*(-(1/2)*R^2/L-(1/2)*(C*(C*R^2-4*L))^(1/2)*R/(C*L))*V__0*(R^2*C^2-C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L+2*(C^2*R^2-4*C*L)^(1/2))*exp(-(1/2)*(R*C+(C*(C*R^2-4*L))^(1/2))*t/(C*L))/(C*R^2-4*L)-(1/2)*(-(1/2)*R^2/L+(1/2)*(C*(C*R^2-4*L))^(1/2)*R/(C*L))*(R^2*C^2+C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L-2*(C^2*R^2-4*C*L)^(1/2))*V__0*exp(-(1/2)*(R*C-(C*(C*R^2-4*L))^(1/2))*t/(C*L))/(C*R^2-4*L), q(t) = -(1/2)*exp(-(1/2)*(R*C-(C*(C*R^2-4*L))^(1/2))*t/(C*L))*(R^2*C^2+C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L-2*(C^2*R^2-4*C*L)^(1/2))*V__0/(C*R^2-4*L)-(1/2)*exp(-(1/2)*(R*C+(C*(C*R^2-4*L))^(1/2))*t/(C*L))*V__0*(R^2*C^2-C*(C^2*R^2-4*C*L)^(1/2)*R-4*C*L+2*(C^2*R^2-4*C*L)^(1/2))/(C*R^2-4*L)+V__0*C}

> 

ansnum := eval(ans, params):
> 

``
> 

plot(eval(V__C(t), ansnum), t = 0 .. 6);

> 

plot(eval(diff(q(t), t), ansnum), t = 0 .. 6);

> 

``


 

Download RLC.mw

Matrix...

dharr 8507
February 06 2022
2 1

If it is a Matrix before you set the values, then you only need to give the name of the matrix to output it,

Matrix.mw

eval, assign...

dharr 8507
February 06 2022
3 0

Use eval for this.

> 

x := 2*c[1]+c[2] = 5;

2*c[1]+c[2] = 5

> 

y := c[1]+c[2] = 3;

c[1]+c[2] = 3

> 

ans := solve({x, y}, {c[1], c[2]})

{c[1] = 2, c[2] = 1}

> 

z := c[1]^2+c[2];

c[1]^2+c[2]

> 

a := eval(z, ans);

5

> 

b := eval(c[1], ans);

2

If you want to set the values of c[1] and c[2] from now on, you can just use assign (but then equations like x and y can't be used as equations again)

> 

assign(ans);
> 

c[1];

2

1

> 

 

Download eval.mw

plottools:-line...

dharr 8507
February 05 2022
0 3

There are a lot of ways to do this. But my advice for the new Maple user would be to use plot3d and other plotting commands from the plots and plottools packages rather than PLOT3D, which is a low-level way for the more technical user. Here is one way:

> 

restart;
> 

with(plots):
n:=2;

2

Unit vectors

> 

a:=[1,0,0]:b:=[0,1,0]:c:=[0,0,1]:

Cube origins

> 

pts:=[seq(seq(seq([i,j,k],i=0..n),j=0..n),k=0..n)]:
 

Cube unit vectors relative to the origins, then remove lines "sticking out"

> 

lines:=map(o->([o,o+a],[o,o+b],[o,o+c]),pts):
lines:=remove(has,lines,n+1):
> 

plotlines:=display(map(x->plottools:-line(x[]),lines),color=blue,linestyle=solid,thickness=2,axes=none):
display(plotlines);

With points as well

> 

display(plotlines,pointplot3d(pts,color=red,symbol=solidsphere,symbolsize=20));

> 

 

Download grid.mw

another way...

dharr 8507
January 27 2022
3 0
> 

restart;
p:=(r*exp(I*theta)-6*I);
eq:=simplify(evalc(p*conjugate(p)))=9; #squared equation

r*exp(I*theta)-6*I

r^2-12*r*sin(theta)+36 = 9

> 

plots:-implicitplot(eq,r=0..10,theta=0..2*Pi,gridrefine=2);

> 

deq:=implicitdiff(eq,theta,r); #want dtheta/dr=0
solve({eq,deq,r>0},{r,theta}); #(could also specify theta<2 to avoid second root)
z:=evalc(eval(r*exp(I*theta),%[1]));

-(1/6)*(-r+6*sin(theta))/(cos(theta)*r)

{r = 3*3^(1/2), theta = (1/3)*Pi}, {r = 3*3^(1/2), theta = (2/3)*Pi}

(3/2)*3^(1/2)+(9/2)*I

> 

 

Download complex.mw

tubeplot...

dharr 8507
January 21 2022
5 1

I don't think you can color a spacecurve with a color function, which is I think what odeplot produces, but a thin tubeplot is near enough. So here is one way.
 

> 

restart;
> 

with(plots):
> 

sys := {
   diff(x(t), t) = v(t)
  ,diff(v(t), t) = cos(t)
  ,x(0) = -1
  ,v(0) = 0

  ,px(t) = piecewise(x(t) >=0, 1, -1)
  ,pv(t) = piecewise(v(t) >=0, 1, -1)
}:
sol := dsolve(sys, [x(t),v(t),px(t),pv(t)], numeric,output=listprocedure): #force variable order with list

Procedures to generate t,x,v,px,pv; take t as argument

> 

T,X,V,PX,PV:=rhs~(sol)[]:

Colours translation table

> 

cols:=table():
cols[1,1]:=[0,255,0]:#"Green":
cols[1,-1]:=[255,0,0]:#"Red":
cols[-1,1]:=[0,0,255]:#"Blue":
cols[-1,-1]:=[255,215,0]:#"Gold":

Colour functions  - takes t and returns red component of required colour; then green then blue

> 

fR := t->cols[round(PX(t)),round(PV(t))][1]:
fG := t->cols[round(PX(t)),round(PV(t))][2]:
fB := t->cols[round(PX(t)),round(PV(t))][3]:
> 

tubeplot([t,X(t),V(t)],t=0..4*Pi,radius=0.05,color=[fR(t),fG(t),fB(t)],style=surface,lightmodel=none,labels=[t,x(t),v(t)]);

> 

 

 

Download coloring2.mw

You can skip px and pv and work directly with x and y for the colours in this case, but I assume you have a more complex case in mind. Version without px or pv:

coloring3.mw

names starting with a digit...

dharr 8507
January 21 2022
1 0

I agree with @nm and @tomleslie that the double underline subscript is better if you just want it to look nice. If you really want an indexed name with index 1d you can enclose it in back quotes: `1d`.

As for why it appears as 1., note that the floating point number 1 can be written as 1e0 or just 1e - see the ?float help page

Surprisingly, 1d does the same as 1e - try 1d3 to get 1000.

I'm guessing this is a holdover from the FORTRAN days, when d would have implied double precision

complexplot3d...

dharr 8507
January 21 2022
2 1

For example

plots:-complexplot3d(Zeta(z),z=-4-10*I..4+40*I,view=[default,default,-1..3]);


Of course you can play around with the ranges to see different parts. This is the magnitude and the colours give the phase.

Download Zeta.mw

version dependent...

dharr 8507
January 16 2022
2 0

Interesting - it worked in version 2015. Probably evalc(abs(...)) would force it.

> 

interface(version)

`Standard Worksheet Interface, Maple 2015.1, Windows 8, June 4 2015 Build ID 1049007`

> 

restart;
> 

v:=1/16*(I*5^(1/2)+I-4*sin(1/5*Pi))*2^(3/10)*(5-5^(1/2))^(1/2)-1/8*2^(4/5)*((I*5^(1/2)+I)*sin(1/5*Pi)+1/2*5^(1/2)+3/2)

(1/16)*(I*5^(1/2)+I-4*sin((1/5)*Pi))*2^(3/10)*(5-5^(1/2))^(1/2)-(1/8)*2^(4/5)*((I*5^(1/2)+I)*sin((1/5)*Pi)+(1/2)*5^(1/2)+3/2)

> 

abs(v)

((-(1/4)*(5-5^(1/2))^(1/2)*2^(3/10)*sin((1/5)*Pi)-(1/8)*2^(4/5)*((1/2)*5^(1/2)+3/2))^2+((1/16)*(5-5^(1/2))^(1/2)*2^(3/10)*(5^(1/2)+1)-(1/8)*2^(4/5)*(5^(1/2)+1)*sin((1/5)*Pi))^2)^(1/2)

> 

sqrt( Re(v)^2 + Im(v)^2 );

((-(1/4)*(5-5^(1/2))^(1/2)*2^(3/10)*sin((1/5)*Pi)-(1/8)*2^(4/5)*((1/2)*5^(1/2)+3/2))^2+((1/16)*(5-5^(1/2))^(1/2)*2^(3/10)*(5^(1/2)+1)-(1/8)*2^(4/5)*(5^(1/2)+1)*sin((1/5)*Pi))^2)^(1/2)

> 

 

Download complex.mw

extended or standard typesetting...

dharr 8507
January 16 2022
1 6

Perhaps this is what you want. The default is standard in earlier versions of Maple, but extended in more recent versions.

Use of Typesetting:-RuleAssistant(); (or view->typesetting from the menu) gives you finer control, or commands from the Typesetting package.

> 

restart;
> 

interface(typesetting=standard):
> 

BesselJ(nu,z);

BesselJ(nu, z)

> 

interface(typesetting=extended):
> 

BesselJ(nu,z);

BesselJ(nu, z)

> 

 

Download typesetting.mw

First 46 47 48 49 50 51 52 Last Page 48 of 83