Advanced Search

dharr

Dr. David Harrington

8235 Reputation

22 Badges

20 years, 340 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

change of variables...

dharr 8235
March 03 2022
7 4

Not sure why, but a change of variables works.
 

> 

restart
> 

ode1 := diff(T__3(t), t, t)+3*(diff(a(t), t))*(diff(T__3(t), t))/a(t)+(2*(diff(a(t), t, t))/a(t)+6*(diff(a(t), t))^2/a(t)^2+(-Omega^2+6)/a(t)^2)*T__3(t)

diff(diff(T__3(t), t), t)+3*(diff(a(t), t))*(diff(T__3(t), t))/a(t)+(2*(diff(diff(a(t), t), t))/a(t)+6*(diff(a(t), t))^2/a(t)^2+(-Omega^2+6)/a(t)^2)*T__3(t)

lower case zeta since Zeta is the zeta function (type it in rather than 2-D pallete that gives Zeta but looks the same)

> 

at := zeta*(1-(1-t/zeta)^2)^(1/2)

zeta*(1-(1-t/zeta)^2)^(1/2)

> 

ode2 := eval(ode1, a(t) = at);

diff(diff(T__3(t), t), t)+3*(1-t/zeta)*(diff(T__3(t), t))/((1-(1-t/zeta)^2)*zeta)+(2*(-(1-t/zeta)^2/((1-(1-t/zeta)^2)^(3/2)*zeta)-1/((1-(1-t/zeta)^2)^(1/2)*zeta))/(zeta*(1-(1-t/zeta)^2)^(1/2))+6*(1-t/zeta)^2/((1-(1-t/zeta)^2)^2*zeta^2)+(-Omega^2+6)/(zeta^2*(1-(1-t/zeta)^2)))*T__3(t)

Change of variables

> 

itr := z = t/zeta-1;

z = t/zeta-1

t = (1+z)*zeta

(diff(diff(T__3(z), z), z))/zeta^2-3*z*(diff(T__3(z), z))/((-z^2+1)*zeta^2)+(2*(-z^2/((-z^2+1)^(3/2)*zeta)-1/((-z^2+1)^(1/2)*zeta))/(zeta*(-z^2+1)^(1/2))+6*z^2/((-z^2+1)^2*zeta^2)+(-Omega^2+6)/(zeta^2*(-z^2+1)))*T__3(z)

> 

sol := dsolve(ode3, T__3(z));

T__3(z) = _C1*LegendreP((-Omega^2+1)^(1/2)-1/2, ((1/2)*I)*15^(1/2), z)/(z^2-1)^(1/4)+_C2*LegendreQ((-Omega^2+1)^(1/2)-1/2, ((1/2)*I)*15^(1/2), z)/(z^2-1)^(1/4)

Change back to t

> 

solt := PDEtools:-dchange(itr, sol, [t], params = {Omega, zeta})

T__3(t) = _C1*LegendreP((-Omega^2+1)^(1/2)-1/2, ((1/2)*I)*15^(1/2), t/zeta-1)/((t/zeta-1)^2-1)^(1/4)+_C2*LegendreQ((-Omega^2+1)^(1/2)-1/2, ((1/2)*I)*15^(1/2), t/zeta-1)/((t/zeta-1)^2-1)^(1/4)

> 

 


 

Download solve3.mw

exponentials...

dharr 8235
March 02 2022
1 0

I'm not sure I understand fully, but the following shows how to get an exponential rise or fall between two values. If you want to upload a worksheet to demonstrate more clearly your problem, use the green up arrow in the editor.
 

> 

restart;

Exponential from max to min or min to max. Assume a form like

> 

y:=a+b*exp(-c*t);

a+b*exp(-c*t)

From max down to min

> 

eqns:=[eval(y,t=0)=ymax,limit(y,t=infinity)=ymin] assuming c::positive;
ans:=solve(eqns,{a,b});
y1:=eval(y,ans);

[a+b = ymax, a = ymin]

{a = ymin, b = -ymin+ymax}

ymin+(-ymin+ymax)*exp(-c*t)

From min up to max

> 

eqns:=[eval(y,t=0)=ymin,limit(y,t=infinity)=ymax] assuming c::positive;
ans:=solve(eqns,{a,b});
y2:=eval(y,ans);

[a+b = ymin, a = ymax]

{a = ymax, b = ymin-ymax}

ymax+(ymin-ymax)*exp(-c*t)

> 

plot(eval([y1,y2],{c=1,ymin=2,ymax=5}),t=0..5,colour=[red,blue]);

> 

 


 

Download exp.mw

plots package...

dharr 8235
February 28 2022
2 1

display is in the plots package, so you need

plots:-display(pp1) ;

(or use with(plots) at the beginning of your worksheet)

odeplot...

dharr 8235
February 28 2022
2 0

odeplot can handle most sorts of plots for numeric odes.

[worksheet display not working right now]

odeplot.mw

Equivalent Resistance...

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

First 44 45 46 47 48 49 50 Last Page 46 of 81