Advanced Search

dharr

Dr. David Harrington

8482 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

Maple-assisted...

dharr 8482
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 8482
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 8482
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 8482
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 8482
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 8482
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 8482
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 8482
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 8482
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 8482
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 8482
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 8482
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

save as eps...

dharr 8482
January 14 2022
2 0

On my windows system and Maple 2015, copying a 2D plot only gives bitmap options under paste special (in CorelDraw). But using the right-click menu to export the plot as an eps file gives a vector graphics file that I can import into CorelDraw and edit. For my version of Maple, this doesn't work with 3D plots, but that might be fixed in later versions.

Generalized eigenvalue problem...

dharr 8482
January 12 2022
3 1

You aren't clear about what you want or what you have tried. If I assume that omega^2 are the eigenvalues (and not the lambdas), then this is a generalized eigenvalue problem for which Maple's Eigenvectors routine gives an answer in a straighforward way. I then forced a change of algorithm whch makes some of the apparently infinite eigenvalues to change to zero, but everything else looks the same.

 

> 

restart;
> 

interface(version);with(LinearAlgebra):

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

> 

M:=4;N:=2;
gen:=generator=-5.0..5.0;
gen2:=generator=1..10; #positive

4

2

generator = -5.0 .. 5.0

generator = 1 .. 10

Generate positive definite matrices

> 

U:=RandomMatrix(M,gen):
v:=RandomVector(M,gen2):
K11:=U.DiagonalMatrix(v).U^+;
IsDefinite(K11,query=positive_definite);

K11 := Matrix(4, 4, {(1, 1) = 474.7658127874699, (1, 2) = -394.5805117919539, (1, 3) = 351.8490779680129, (1, 4) = 198.4791447697421, (2, 1) = -394.5805117919539, (2, 2) = 360.93618001079045, (2, 3) = -217.9024802412394, (2, 4) = -174.93938764201536, (3, 1) = 351.8490779680129, (3, 2) = -217.9024802412394, (3, 3) = 505.123761706797, (3, 4) = 265.84210885475943, (4, 1) = 198.4791447697421, (4, 2) = -174.93938764201536, (4, 3) = 265.84210885475943, (4, 4) = 470.4006207069442}, datatype = float[8])

true

> 

U:=RandomMatrix(M,gen):
v:=RandomVector(M,gen2):
M11:=U.DiagonalMatrix(v).U^+;
IsDefinite(M11,query=positive_definite);

M11 := Matrix(4, 4, {(1, 1) = 272.4454727874966, (1, 2) = 22.89942145896406, (1, 3) = 159.03066901418603, (1, 4) = -83.8538688915777, (2, 1) = 22.899421458964053, (2, 2) = 68.23916321060702, (2, 3) = 99.66785764378282, (2, 4) = 53.80202849499174, (3, 1) = 159.030669014186, (3, 2) = 99.66785764378282, (3, 3) = 274.1189955623002, (3, 4) = 48.570773916747285, (4, 1) = -83.8538688915777, (4, 2) = 53.80202849499174, (4, 3) = 48.57077391674727, (4, 4) = 158.7300008588166}, datatype = float[8])

true

Assemble K and M matrices

> 

K12:=RandomMatrix(M,N,gen):
KK1:=Matrix(<<K11|K12>,<K12^+|Matrix(N,N,0)>>,shape=symmetric);
IsDefinite(KK1,query=positive_semidefinite);

KK1 := Matrix(6, 6, {(1, 1) = HFloat(474.7658127874699), (1, 2) = HFloat(-394.5805117919539), (1, 3) = HFloat(351.8490779680129), (1, 4) = HFloat(198.4791447697421), (1, 5) = 2.95199901137063137, (1, 6) = -4.65553919497091240, (2, 2) = HFloat(360.93618001079045), (2, 3) = HFloat(-217.9024802412394), (2, 4) = HFloat(-174.93938764201536), (2, 5) = 2.65516788149002370, (2, 6) = 4.50222048838354993, (3, 3) = HFloat(505.123761706797), (3, 4) = HFloat(265.84210885475943), (3, 5) = -1.18441542906991604, (3, 6) = -1.82900519939139450, (4, 4) = HFloat(470.4006207069442), (4, 5) = -.612556403436017582, (4, 6) = 1.94828622975817023, (5, 5) = 0, (5, 6) = 0, (6, 6) = 0}, storage = triangular[upper], shape = [symmetric])

false

> 

MM1:=Matrix(<<M11|Matrix(M,N,0)>,<Matrix(N,M,0)|Matrix(N,N,0)>>,shape=symmetric);
IsDefinite(MM1,query=positive_semidefinite);

MM1 := Matrix(6, 6, {(1, 1) = HFloat(272.4454727874966), (1, 2) = HFloat(22.89942145896406), (1, 3) = HFloat(159.03066901418603), (1, 4) = HFloat(-83.8538688915777), (1, 5) = 0, (1, 6) = 0, (2, 2) = HFloat(68.23916321060702), (2, 3) = HFloat(99.66785764378282), (2, 4) = HFloat(53.80202849499174), (2, 5) = 0, (2, 6) = 0, (3, 3) = HFloat(274.1189955623002), (3, 4) = HFloat(48.570773916747285), (3, 5) = 0, (3, 6) = 0, (4, 4) = HFloat(158.7300008588166), (4, 5) = 0, (4, 6) = 0, (5, 5) = 0, (5, 6) = 0, (6, 6) = 0}, storage = triangular[upper], shape = [symmetric])

true

Solve the generalized eigenvector problem. Shape=symmetric is not enough to force a real calculation and give floating eigenvalues, but get complex eigenvalueswith zero imaginary parts.

> 

a,A:=Eigenvectors(KK1,MM1):
a:=simplify(a,zero); #omega^2
A:=simplify(A,zero); #Matrix of eigenvectors

a := Vector(6, {(1) = Float(infinity), (2) = 6.793175977, (3) = .9218354408, (4) = Float(infinity), (5) = Float(-infinity), (6) = Float(infinity)})

A := Matrix(6, 6, {(1, 1) = -0.1839610129e-16, (1, 2) = 0.2525244976e-2, (1, 3) = 0.2283298399e-2, (1, 4) = -0.1082838565e-14, (1, 5) = 0.1046046362e-14, (1, 6) = 0.1082838565e-14, (2, 1) = -0.1268034289e-16, (2, 2) = -0.1710091008e-2, (2, 3) = -0.1599152764e-1, (2, 4) = 0.2921824402e-15, (2, 5) = -0.3175431260e-15, (2, 6) = -0.2921824402e-15, (3, 1) = 0.7177805376e-17, (3, 2) = -0.1820461508e-2, (3, 3) = -0.3506659550e-1, (3, 4) = 0.5808392185e-15, (3, 5) = -0.5664836078e-15, (3, 6) = -0.5808392185e-15, (4, 1) = 0.3540019806e-18, (4, 2) = 0.8276992328e-2, (4, 3) = 0.9490588590e-2, (4, 4) = -0.5375652816e-15, (4, 5) = 0.5382732856e-15, (4, 6) = 0.5375652816e-15, (5, 1) = -1., (5, 2) = -.1674559135, (5, 3) = -1., (5, 4) = -1., (5, 5) = -1., (5, 6) = 1., (6, 1) = .8473012185, (6, 2) = 1., (6, 3) = -0.8033014586e-1, (6, 4) = .8473012185, (6, 5) = .8473012185, (6, 6) = -.8473012185})

Force different algorithm by declaring M with positive definite attribute (see Matrix help page). But since it isn't this might not work.
Find the eigenvalues are now real, and some of the infinite ones are now zero. Eigenvectors are just the same.

> 

MM2:=Matrix(MM1,shape=symmetric,attributes=[positive_definite]):
> 

a2,A2:=Eigenvectors(KK1,MM2):
a2;
simplify(A2,zero);

Vector[column]([[HFloat(HFloat(infinity))], [0.], [6.79317597683730], [0.], [.921835440824920], [0.]])

Matrix([[-0.1411084663e-17, 0.2525244976e-2, 0.2283298399e-2, -0.1065853548e-14, 0.1063031379e-14, 0.1065853548e-14], [-0.6950197323e-18, -0.1710091008e-2, -0.1599152764e-1, 0.3041677633e-15, -0.3055578028e-15, -0.3041677633e-15], [0.3934212517e-18, -0.1820461508e-2, -0.3506659550e-1, 0.5740548344e-15, -0.5732679919e-15, -0.5740548344e-15], [0.1940313160e-19, 0.8276992328e-2, 0.9490588590e-2, -0.5378998805e-15, 0.5379386867e-15, 0.5378998805e-15], [-1., -.1674559135, -1., -1., -1., 1.], [.8473012185, 1., -0.8033014586e-1, .8473012185, .8473012185, -.8473012185]])

> 

 

 

Download evecs4.mw

Edit: noticed the eigenvalues change order but the eigenvectors did not. Interesting.

Eigenvectors needs shape=symmetric for both matrices and only the positive_definite attribute for the second one to produce real eigevalues via the faster algorithm.

evalc...

dharr 8482
January 09 2022
6 0

As @tomleslie says, uploading a worksheet would be helpful. I'm guessing Maple is lacking some information to proceed further, or perhaps you used a=expression rather than a:=expression. evalc can give a value for abs assuming the variables are real:

> 

a := 1/(1+I*omega*tau)

1/(1+I*omega*tau)

Maple doesn't know enough about omega and tau to simplify any further than

> 

abs(a);

1/abs(1+I*omega*tau)

evalc forces evaluation assuming all the indeterminates are real

> 

evalc(abs(a));

1/(omega^2*tau^2+1)^(1/2)

Or the assumptions can be explicitly stated

> 

`assuming`([abs(a)], [positive]);

1/(omega^2*tau^2+1)^(1/2)

> 

 

Download evalc.mw

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