dharr

Dr. David Harrington

8205 Reputation

22 Badges

20 years, 337 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
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


 

restart; #dharr comments in red

# Assuming 1/3<qc,h<1 (h does not appear below)and deltac+deltah=1, how to write "w" as a fraction like the fraction below:
# w=(M^2-M1^2)/(M^2*(M^2-M0^2)) where M1 and M0 are functions of (qc,qh,Tch,alpha,deltac,deltah)
#
V__phi:=-(1-((M**2)/((deltac)*(1+(qc-1)*phi)**((qc+1)/(2*qc-2))+(deltah)*(1+(qh-1)*phi)**((qh+1)/(2*qh-2)))**3)*((deltac/2)*(qc+1)*(1+(qc-1)*phi)**((3-qc)/(2*qc-2))+(deltah/2)*Tch*(qh+1)*(1+(qh-1)*phi)**((3-qh)/(2*qh-2))))**(-2)*(((1+(alpha/M)**2*(2*deltac/(3*qc-1)+2*deltah/(Tch*(3*qh-1)))))*(2*deltac/(3*qc-1)*(1+(qc-1)*phi)**((3*qc-1)/(2*qc-2))+2*deltah/(Tch*(3*qh-1))*(1+Tch*(qh-1)*phi)**((3*qh-1)/(2*qh-2))-2*deltac/(3*qc-1)-2*deltah/(Tch*(3*qh-1)))-phi-(1/2)*(alpha/M)**2*((2*deltac/(3*qc-1)*(1+(qc-1)*phi)**((3*qc-1)/(2*qc-2))+2*deltah/(Tch*(3*qh-1))*(1+Tch*(qh-1)*phi)**((3*qh-1)/(2*qh-2)))**2-(2*deltac/(3*qc-1)+2*deltah/(Tch*(3*qh-1)))**2)+(M**2*((1+(alpha/M)**2*(2*deltac/(3*qc-1)+2*deltah/(Tch*(3*qh-1))))))*((((deltac*(1+(qc-1)*phi)**((qc+1)/(2*qc-2)))+(deltah*(1+Tch*(qh-1)*phi)**((qh+1)/(2*qh-2)))))**(-1)-1/(deltac+deltah))-(M**2/2)*((((deltac*(1+(qc-1)*phi)**((qc+1)/(2*qc-2)))+(deltah*(1+Tch*(qh-1)*phi)**((qh+1)/(2*qh-2)))))**(-2)-(deltac+deltah)**(-2))+alpha**2*phi-alpha**2*(((2*deltac/(3*qc-1)*(1+(qc-1)*phi)**((3*qc-1)/(2*qc-2))+2*deltah/(Tch*(3*qh-1))*(1+Tch*(qh-1)*phi)**((3*qh-1)/(2*qh-2))))/((((deltac*(1+(qc-1)*phi)**((qc+1)/(2*qc-2)))+(deltah*(1+Tch*(qh-1)*phi)**((qh+1)/(2*qh-2))))))-((2*deltac/(3*qc-1)+2*deltah/(Tch*(3*qh-1))))/(deltac+deltah))):
#
# use deltac+deltah=1 to eliminate deltah and simplify
w:=simplify(eval(diff(V__phi,phi$2),{phi=0,deltah=1-deltac}));

((((-qh-1)*Tch+qc+1)*deltac+Tch*(qh+1))*M^2-2*alpha^2)/(M^2*(-2+(((-qh-1)*Tch+qc+1)*deltac+Tch*(qh+1))*M^2))

w2:=(M^2-M1^2)/(M^2*(M^2-M0^2));

(M^2-M1^2)/(M^2*(M^2-M0^2))

A few possibilities that might be narrowed down by working with the assumptions.

solve(identity(w=w2,M),{M0,M1});

{M0 = -2/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2), M1 = -2*alpha/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2)}, {M0 = 2/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2), M1 = 2*alpha/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2)}, {M0 = -2/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2), M1 = 2*alpha/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2)}, {M0 = 2/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2), M1 = -2*alpha/(-2*Tch*deltac*qh-2*Tch*deltac+2*Tch*qh+2*deltac*qc+2*Tch+2*deltac)^(1/2)}

NULL

Download A1.mw

Just use Profile(PickAngles);
Then use PickAngles (one or more times)

Then use PrintProfiles(PickAngles)

An example is given here (in the last item of the thread)

The usual way to get math output is to use a Mathematical Expression component (also called a Math Container) and set the expression property:

DocumentTools:-SetProperty("MathContainer0", expression, ifactor(2600))

MathExpression.mw

I don't know much about this, but I understand from the wikipedia page that we just multiply the series by the appropriate root of unity. The following code is not very robust because it looks only at the first term of the series that puiseux outputs. If this is what you want, then it could be improved.

Edit: turns out each term gets a different phase factor, corrected in later post.

restart

with(algcurves); f := 6*w^14+z^30+z^32+w^15*(z^2+2)+w^12*(z^3+z)+w^9*(z^9+z^5)+w^5*(z^20+z^14); puiseuxList := convert(puiseux(f, z = 0, w, 5), list)

6*w^14+z^30+z^32+w^15*(z^2+2)+w^12*(z^3+z)+w^9*(z^9+z^5)+w^5*(z^20+z^14)

[(-z)^(9/4), -(-z)^(16/5), -3-(350887/472392)*z^4-(4381/52488)*z^3+(365/243)*z^2-(1/18)*z, -16*(-z)^(14/3)+(2/3)*(-z)^(13/3)+(1/3)*(-z)^(10/3)+2*z^3-(-z)^(4/3), (35663903797/2166612408926208)*(-6*z)^(9/2)-(3407/944784)*z^4-(28431487/69657034752)*(-6*z)^(7/2)-(310547/104976)*z^3-(62285/26873856)*(-6*z)^(5/2)-(1/972)*z^2-(5/15552)*(-6*z)^(3/2)+(1/36)*z-(1/6)*(-6*z)^(1/2)]

p := puiseuxList[4]

-16*(-z)^(14/3)+(2/3)*(-z)^(13/3)+(1/3)*(-z)^(10/3)+2*z^3-(-z)^(4/3)

We assume the ramification index is the denominator of the exponent of the first term

conj:=proc(p,z,i)
 local term1:=op(1,p),
       ram:=denom(diff(term1,z)/term1*z);
 p*exp(i*2*Pi*I/ram);
end proc;

proc (p, z, i) local term1, ram; term1 := op(1, p); ram := denom((diff(term1, z))*z/term1); p*exp((2*I)*i*Pi/ram) end proc

p1 := conj(p, z, 1); p2 := conj(p, z, 2)

(-16*(-z)^(14/3)+(2/3)*(-z)^(13/3)+(1/3)*(-z)^(10/3)+2*z^3-(-z)^(4/3))*(-1/2+((1/2)*I)*3^(1/2))

(-16*(-z)^(14/3)+(2/3)*(-z)^(13/3)+(1/3)*(-z)^(10/3)+2*z^3-(-z)^(4/3))*(-1/2-((1/2)*I)*3^(1/2))

pa := plot3d([Re(r*exp(I*t)), Im(r*exp(I*t)), Im(eval(p, z = r*exp(I*t)))], r = 0 .. .15, t = 0 .. 2*Pi)

pb := plot3d([Re(r*exp(I*t)), Im(r*exp(I*t)), Im(eval(p1, z = r*exp(I*t)))], r = 0 .. .15, t = 0 .. 2*Pi)

pc := plot3d([Re(r*exp(I*t)), Im(r*exp(I*t)), Im(eval(p2, z = r*exp(I*t)))], r = 0 .. .15, t = 0 .. 2*Pi)

NULL

Download testWorkSheet.mw

The easiest way to get material onto this site is to upload a whole worksheet file using the green up-arrow. Other comments below.

restart

h := proc (z) options operator, arrow; z^(1/2) end proc

proc (z) options operator, arrow; z^(1/2) end proc

plot3d(Im(h(x+I*y)), x = -1 .. 1, y = -1 .. 1)

Changed to a list since order matters. (= changed to :=)

functionList := [2+z, z^2-3*z, -z^3+4]

[2+z, z^2-3*z, -z^3+4]

use unapply to turn an existing expression into a function since the z in v := z -> functionList[1]; is local to v and not the same as the z in the functionList

(I made a 2 argument fuction so the function can also be selected

v := unapply(functionList[i], i, z)

proc (i, z) options operator, arrow; [z+2, z^2-3*z, -z^3+4][i] end proc

v(1, z); v(2, z)

2+z

z^2-3*z

plot3d(Im(v(1, x+I*y)), x = -1 .. 1, y = -1 .. 1)

NULL

Download unappply.mw

So you want a random one of the permutations of 3 letters, which can be done without generating them all as follows:

restart;

it := Iterator:-Permute(26, 3); # permutations of 3 things out of 26
num := Number(it); # number of permutations
rank := rand(1 .. num)(); # randomly choose one
perm:=Unrank(it,rank); # find out what it is
A11, A12, A13 := map(convert, [seq(StringTools:-ExpandCharacterClass("a-z")[i], i in perm)], name)[];

Permute(26, 3)

num := 15600

rank := 3003

Array(%id = 36893490298184220908)

f, a, d

NULL

Download random.mw

I did it in exponential form, looking for exact roots and then finding numerical values for them.

Edit: I redid it in cartesian form. Earlier exponential form below

Your error message for RootFinding:-Isolate symbolic coefficients are not supported, {f, i, k} means you need to give values to f,i,and k or use solve for symbolic solutions.

restart;

Solve the following for x,y,z complex of magnitude 1 and z a d-root of unity

p:=y*(1-x^(m+1)*z)+(1-x^(n+1)*z);

y*(1-x^(m+1)*z)+1-x^(n+1)*z

Take the case of m=1, n=3 and d=3 or z^3=1

p133:=eval(p,{m=1,n=3});

y*(-x^2*z+1)+1-x^4*z

Choose primitive root for z for example. Set up 4 eqns in 4 unknowns

k:=1; # or k=0 or k=2
xyz[k]:={x=a+I*b,y=c+d*I,z=evalc(exp(k*2*Pi*I/3))};
eq:=evalc(eval(p133,xyz[k])):
eqRe:=evalc(Re(eq));eqIm:=evalc(Im(eq));
eqx:=a^2+b^2=1;eqy:=c^2+d^2=1;

1

{x = a+I*b, y = c+I*d, z = -1/2+((1/2)*I)*3^(1/2)}

1+2*3^(1/2)*(a^3*b-a*b^3)+c*((1/2)*a^2-(1/2)*b^2+a*b*3^(1/2)+1)-d*(a*b+(1/2)*(-a^2+b^2)*3^(1/2))+(1/2)*a^4-3*a^2*b^2+(1/2)*b^4

d*((1/2)*a^2-(1/2)*b^2+a*b*3^(1/2)+1)+c*(a*b+(1/2)*(-a^2+b^2)*3^(1/2))+2*a^3*b-2*a*b^3+(1/2)*(-a^4+6*a^2*b^2-b^4)*3^(1/2)

a^2+b^2 = 1

c^2+d^2 = 1

Find 8 solutions

sols[k]:=simplify([solve({eqRe,eqIm,eqx,eqy},{a,b,c,d},explicit)]);
nops(sols[k]);

[{a = 1, b = 0, c = -1, d = 0}, {a = -1, b = 0, c = -1, d = 0}, {a = (1/4)*((4+(4*I)*3^(1/2))^(2/3)+4)/(4+(4*I)*3^(1/2))^(1/3), b = (I*(4+(4*I)*3^(1/2))^(1/3)-I+3^(1/2))/(4+(4*I)*3^(1/2))^(2/3), c = -(1/4)*((4+(4*I)*3^(1/2))^(2/3)+4)/(4+(4*I)*3^(1/2))^(1/3), d = -(I*(4+(4*I)*3^(1/2))^(1/3)-I+3^(1/2))/(4+(4*I)*3^(1/2))^(2/3)}, {a = -3^(1/2)*sin((4/9)*Pi)+cos((1/9)*Pi), b = -(1/2)*((I-3^(1/2))*(4+(4*I)*3^(1/2))^(1/3)-4*I)/(4+(4*I)*3^(1/2))^(2/3), c = 3^(1/2)*sin((4/9)*Pi)-cos((1/9)*Pi), d = (1/2)*((I-3^(1/2))*(4+(4*I)*3^(1/2))^(1/3)-4*I)/(4+(4*I)*3^(1/2))^(2/3)}, {a = 3^(1/2)*sin((4/9)*Pi)-2*cos((1/9)*Pi), b = -(1/2)*((4+(4*I)*3^(1/2))^(1/3)+2)*(3^(1/2)+I)/(4+(4*I)*3^(1/2))^(2/3), c = -3^(1/2)*sin((4/9)*Pi)+2*cos((1/9)*Pi), d = (1/2)*((4+(4*I)*3^(1/2))^(1/3)+2)*(3^(1/2)+I)/(4+(4*I)*3^(1/2))^(2/3)}, {a = (1/4)*((-4+(4*I)*3^(1/2))^(2/3)+4)/(-4+(4*I)*3^(1/2))^(1/3), b = -(I*(-4+(4*I)*3^(1/2))^(1/3)+3^(1/2)+I)/(-4+(4*I)*3^(1/2))^(2/3), c = (1/4)*((-4+(4*I)*3^(1/2))^(2/3)+4)/(-4+(4*I)*3^(1/2))^(1/3), d = -(I*(-4+(4*I)*3^(1/2))^(1/3)+3^(1/2)+I)/(-4+(4*I)*3^(1/2))^(2/3)}, {a = -cos((1/9)*Pi), b = (1/2)*((-4+(4*I)*3^(1/2))^(1/3)-2)*(I-3^(1/2))/(-4+(4*I)*3^(1/2))^(2/3), c = -cos((1/9)*Pi), d = (1/2)*((-4+(4*I)*3^(1/2))^(1/3)-2)*(I-3^(1/2))/(-4+(4*I)*3^(1/2))^(2/3)}, {a = -3^(1/2)*sin((4/9)*Pi)+2*cos((1/9)*Pi), b = (1/2)*((3^(1/2)+I)*(-4+(4*I)*3^(1/2))^(1/3)+4*I)/(-4+(4*I)*3^(1/2))^(2/3), c = -3^(1/2)*sin((4/9)*Pi)+2*cos((1/9)*Pi), d = (1/2)*((3^(1/2)+I)*(-4+(4*I)*3^(1/2))^(1/3)+4*I)/(-4+(4*I)*3^(1/2))^(2/3)}]

8

Convert to 8 solutions for x,y

solsxy[k]:=simplify(map2(eval,xyz[k],sols[k]));

[{x = 1, y = -1, z = -1/2+((1/2)*I)*3^(1/2)}, {x = -1, y = -1, z = -1/2+((1/2)*I)*3^(1/2)}, {x = 2*(1+I*3^(1/2))/(4+(4*I)*3^(1/2))^(2/3), y = -2*(1+I*3^(1/2))/(4+(4*I)*3^(1/2))^(2/3), z = -1/2+((1/2)*I)*3^(1/2)}, {x = -4/(4+(4*I)*3^(1/2))^(2/3), y = 4/(4+(4*I)*3^(1/2))^(2/3), z = -1/2+((1/2)*I)*3^(1/2)}, {x = -2*(-1+I*3^(1/2))/(4+(4*I)*3^(1/2))^(2/3), y = 2*(-1+I*3^(1/2))/(4+(4*I)*3^(1/2))^(2/3), z = -1/2+((1/2)*I)*3^(1/2)}, {x = 2/(-4+(4*I)*3^(1/2))^(1/3), y = 2/(-4+(4*I)*3^(1/2))^(1/3), z = -1/2+((1/2)*I)*3^(1/2)}, {x = -(1/2)*(4+(4*I)*3^(1/2))^(1/3), y = -(1/2)*(4+(4*I)*3^(1/2))^(1/3), z = -1/2+((1/2)*I)*3^(1/2)}, {x = (4+(4*I)*3^(1/2))^(1/3)*(-1+I*3^(1/2))/((2*I)*3^(1/2)+2), y = (4+(4*I)*3^(1/2))^(1/3)*(-1+I*3^(1/2))/((2*I)*3^(1/2)+2), z = -1/2+((1/2)*I)*3^(1/2)}]

Check

simplify(map2(eval,p133,solsxy[k]));

[0, 0, 0, 0, 0, 0, 0, 0]

Numerical values

evalf(solsxy[k]);

[{x = 1., y = -1., z = -.5000000000+.8660254040*I}, {x = -1., y = -1., z = -.5000000000+.8660254040*I}, {x = .9396926206+.3420201432*I, y = -.9396926206-.3420201432*I, z = -.5000000000+.8660254040*I}, {x = -.7660444428+.6427876096*I, y = .7660444428-.6427876096*I, z = -.5000000000+.8660254040*I}, {x = -.1736481779-.9848077528*I, y = .1736481779+.9848077528*I, z = -.5000000000+.8660254040*I}, {x = .7660444432-.6427876096*I, y = .7660444432-.6427876096*I, z = -.5000000000+.8660254040*I}, {x = -.9396926210-.3420201433*I, y = -.9396926210-.3420201433*I, z = -.5000000000+.8660254040*I}, {x = .1736481779+.9848077531*I, y = .1736481779+.9848077531*I, z = -.5000000000+.8660254040*I}]

NULL

Download solvecartesian.mw

[Exponential form: Download solve.mw]

The proposed solution satisfies the PDE:

Edit: you had some strange 2-D code that wasn't working correctly (no sqrt) so I originally thought it was not a solution.

solution_pde_check.mw

In general, plots:-display can combine various plots,so it can be done this way.

pde.mw

You write your application in Maple using embedded components that you interact with only by usng those components - say by pushing the button. Your plot would be produced in a plot component. Then you just open it in Maple Player and it works in exactly the same way.

You do not need to upgrade Maple to use the latest Maple Player. Maple Player is free.

But I don't think it is available for the ipad; the documentation says Mac, Windows and Linux.

I changed pi to Pi (pi has no special meaning in Maple), and sin(2*x) to sin(2*x[i])

The whole x array (each entry multiplied by 2) was being passed to sin.

 Asst_2_Q3bcd.mw

(For the test matrix below, the circles have radius epsilon, as is true for any normal matrix.)

Pseudsospectrum of a matrix - see Wikipedia - Pseudospectrum

restart

A := LinearAlgebra:-DiagonalMatrix([-1., 1.+1.*I, -I])

Matrix(%id = 36893490936251856644)

Contours are at the selected epsilon values; crosses at the eigenvalues

pseudoplot:=proc(A::Matrix(square),epsvals::list(positive))
  local evs,conts,n,maxc,rng,evplot,f;
  uses LinearAlgebra;
  n:=Dimension(A)[1];
  evs:=map([Re,Im],Eigenvalues(evalf(A),'output'='list'));
  conts:=sort(epsvals);
  maxc:=conts[-1];
  rng:=max(map(abs,ListTools:-FlattenOnce(evs)))+maxc;
  evplot:=plots:-pointplot(evs,symbol='diagonalcross');
  f:=(x,y)->1/min(SingularValues((x+I*y)*IdentityMatrix(n)-A));
  plots:-display(
     plots:-contourplot(f, -rng .. rng, -rng .. rng, 'contours' = conts^~(-1)),
     evplot,
     'labels'=["Re","Im"],'axes'='boxed','font'=['times','roman',13],'labelfont'=['times','roman',14]);
end proc:

pseudoplot(A, [2, .2, .6, 1.4])

NULL

Download pseudospectrum.mw

Here's one way to do it. I'm assuming a list of lists is OK.

(Worksheet not displaying right now.)

Download powers.mw

powers:=proc(p,x,y)
 local pow:=(term,x,y)->[diff(term,x)/term*x,diff(term,y)/term*y];
 if type(p,`+`) then
   map(pow,[op(p)],x,y)
 else
   [pow(p,x,y)]
 end if;
end proc:

 

Here's one way to do it.

pointplot.mw

The loop can be done more elegantly generating the sequence (in 1-D) with a for loop

pts := [(for i from 0.1 to 0.9 by 0.1 do
         val := FirmModelPartial1(i, 0.2, 0.1)[3];
         if type(val,numeric) then [i, val] else NULL end if end do)];

but I can't remember if this was possible or not in Maple 2021.

Maple got stuck on the solution with arbitrary constants for longer than I was willing to wait. The equation is separable, and I'm guessing evaluation of the integral needs some appropriate assumptions about the values of the constants. So you can make progress for some specific values of the constants, for example if they are all 1:

restart;

AAS := C*diff(y(t), t) + (-B0*y(t)^3 - B1*y(t)^2 - B2*y(t) - B3);
IC:=y(0)= 1;

C*(diff(y(t), t))-B0*y(t)^3-B1*y(t)^2-B2*y(t)-B3

y(0) = 1

AASnum:=eval(AAS,{C=1,B0=1,B1=1,B2=1,B3=1});

diff(y(t), t)-y(t)^3-y(t)^2-y(t)-1

yt:=dsolve({AASnum,IC});

y(t) = RootOf(8*t+2*ln(_Z^2+1)-4*arctan(_Z)-4*ln(_Z+1)+2*ln(2)+Pi)

plot(rhs(yt),t=0..0.218);

Or, equivalently

sol:=dsolve({AASnum,IC},implicit);

t+(1/4)*ln(y(t)^2+1)-(1/2)*arctan(y(t))-(1/2)*ln(y(t)+1)+(1/4)*ln(2)+(1/8)*Pi = 0

plots:-implicitplot(subs(y(t)=y,sol),t=0..5,y=0..100);

NULL

Download dsolve.mw

First 26 27 28 29 30 31 32 Last Page 28 of 81