Advanced Search

dharr

Dr. David Harrington

8270 Reputation

22 Badges

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

some things...

dharr 8270
October 14 2023
1 8

[Edit: the OP has now completely changed this post, so my responses here and below are no longer relevant.]

 

The equations in your post outline a laplace transform analytical solution and you seem to want to follow that, but in your worksheet you set up for a numerical solution of the original untransformed equations, abandoning the laplace transform method.

So forgetting laplace transforms and following your worksheet, you have two errors in formulating a consistent problem. You had strange parentheses after alpha0 in your Cond, which I guessed as a correction (in red). [Edit] What about initial conditions? 

new_paper_dust.mw

If you want the analytical solution then you can start by transforming the pdes to odes with intrans:-laplace, gving you Eq (16) etc, then solve the ode system - but then it is not a good idea to give numerical values to your parameters.

type of file...

dharr 8270
October 12 2023
1 7

Maple tries to figure out the file format you want from the file extension, but doesn't recognize .dat. You can use the format option to specify which format you want, for example CSV:

> 

restart
> 

vv := Vector([1, 3, 5, 7])

Vector[column](%id = 36893490061895775284)

File output has the four numbers, one per line

> 

Export("c:/Users/dharr/Desktop/testfile.dat", vv, format = "CSV")

8

> 

NULL

Download testfile.mw

Edit: for saving Maple objects in a .m file for reading in later into Maple use the save command

hard to beat...

dharr 8270
October 12 2023
2 3

@sursumCorda Vote up. I could not beat your code, which as you no doubt realize is very efficient because (1) CharacterFrequencies uses external code, (2) it works directly on the strings and (3) anyway most of the time is taken working out the primes.

I tried converting the integer to an Array of digits, and then invoking a compiled routine to search the Array. But the additional conversion time offsets potential efficiency gain, and the prime generation time still dominates. (Actually using numboccur 10 times, dopping out once you find 4, is more efficient and is close to your result)

Given that manipulating digits in integers is probably a common operation of interest I'm surprised that there isn't a better way to extract them. convert(p,base,10) is shockingly slow compared to the methods that go via sprintf("%d",p). Of course the sprint code could probably be easily modified to produce an Array or List.

I'm confused about the relative merits of nextprime over multiple calls to ithprime, and there my timings weren't consistent - perhaps there is an internal list that is added to as another prime is found, or there is a huge table anyway. 

> 

restart;

Compilable procedure to find if any elements in Array e appear exactly four times.

> 

has4:=proc(n::integer,e::Array(datatype=integer[4]),c::Array(datatype=integer[4]))
  local i,ei;
  for i from 1 to 10 do
    c[i]:=0;
  end do;
  for i to n do;
    ei:=e[i]-47;
    c[ei]:=c[ei]+1;
  end do;
  for i from 1 to 10 do
    if c[i]=4 then return true end if;
  end do;
  false;
end proc:
> 

has4c:=Compiler:-Compile(has4);

proc () options call_external, define_external(maple_compiled_m7dc41b208b0ac686076b31608b0d668e, MAPLE, LIB = "C:\Users\dharr\AppData\Local\Temp\dharr-3388\maple_compiled_m7dc41b208b0ac686076b31608b0d668eyYxVbxdp.dll"); call_external(0, 140730599545072, true, false, false, args) end proc

> 

Aseq:=proc(m::nonnegint,n::posint:=1);
  local p,e,ne,p4:=Array([]),c:=Array(1..10,datatype=integer[4]);
  p:=ithprime(n);
  to m do
    e:=Array(convert(sprintf("%d",p),'bytes'),datatype=integer[4]);
    ne:=upperbound(e);
    if has4c(ne,e,c) then p4 ,= p end if;
    p:=nextprime(p);
  end do;
  p4
end proc:
> 

a:=CodeTools:-Usage(Aseq(10^6,10^6)):numelems(a);

memory used=3.40GiB, alloc change=12.00MiB, cpu time=39.20s, real time=34.18s, gc time=7.97s

42720

> 

CodeTools:-Profiling:-Profile(Aseq);
numelems(Aseq(10^6,10^6));

42720

> 

CodeTools:-Profiling:-PrintProfiles(Aseq);

Aseq
Aseq := proc(m::nonnegint, n::posint := 1)
local p, e, ne, p4, c;
     |Calls Seconds  Words|
PROC |    1  44.625 454388814|
   1 |    1   0.000     18| p4 := Array([]);
   2 |    1   0.000     24| c := Array(1 .. 10,datatype = integer[4]);
   3 |    1   0.000     10| p := ithprime(n);
   4 |    1   3.247      7| to m do
   5 |1000000   7.347 55008673|     e := Array(convert(sprintf("%d",p),'bytes'),
                                  datatype = integer[4]);
   6 |1000000   1.490 2000872|     ne := upperbound(e);
   7 |1000000   2.005 4000000|     if has4c(ne,e,c) then
   8 |42720   0.030 127264|         p4 ,= p
                                end if;
   9 |1000000  30.506 393251946|     p := nextprime(p)
                            end do;
  10 |    1   0.000      0| p4
end proc
 

Simpler and unexpectedly faster than the above

> 

restart;
> 

Aseq2:=proc(m::nonnegint,n::posint:=1);
  local p,e,p4:=Array([]);
  p:=ithprime(n);
  to m do
    e:=Array(convert(sprintf("%d",p),'bytes'));
    if numboccur(e,48)=4 or numboccur(e,49)=4 or numboccur(e,50)=4 or numboccur(e,51)=4 or numboccur(e,52)=4 or
       numboccur(e,53)=4 or numboccur(e,54)=4 or numboccur(e,55)=4 or numboccur(e,56)=4 or numboccur(e,57)=4
    then
      p4 ,= p
    end if;
    p:=nextprime(p);
  end do;
  p4
end proc:
> 

a:=CodeTools:-Usage(Aseq2(10^6,10^6)):numelems(a);

memory used=3.80GiB, alloc change=2.00MiB, cpu time=30.56s, real time=30.58s, gc time=2.48s

42720

> 

CodeTools:-Profiling:-Profile(Aseq2);
numelems(Aseq2(10^6,10^6));

42720

> 

CodeTools:-Profiling:-PrintProfiles(Aseq2);

Aseq2
Aseq2 := proc(m::nonnegint, n::posint := 1)
local p, e, p4;
     |Calls Seconds  Words|
PROC |    1  50.343 511132410|
   1 |    1   0.000     18| p4 := Array([]);
   2 |    1   0.000     10| p := ithprime(n);
   3 |    1   2.616      7| to m do
   4 |1000000   7.437 59084011|     e := Array(convert(sprintf("%d",p),'bytes'));
   5 |1000000   6.116 58645626|     if numboccur(e,48) = 4 or numboccur(e,49)
                                  = 4 or numboccur(e,50) = 4 or
                                  numboccur(e,51) = 4 or numboccur(e,
                                  52) = 4 or numboccur(e,53) = 4 or
                                  numboccur(e,54) = 4 or numboccur(e,
                                  55) = 4 or numboccur(e,56) = 4 or
                                  numboccur(e,57) = 4 then
   6 |42720   0.030 127264|         p4 ,= p
                                end if;
   7 |1000000  34.144 393275474|     p := nextprime(p)
                            end do;
   8 |    1   0.000      0| p4
end proc
 

@sursumCorda's routine

> 

restart;
> 

A161786__1 := proc(m::nonnegint, n::posint := 1, ` $`)::'Vector'(prime):
        options no_options:
        local p := ifelse(n = 1, 1, ithprime(n - 1)), vec := DEQueue();
        to m do
                if member(4, rhs~({StringTools['CharacterFrequencies'](nprintf("%d", (p := nextprime(p))), 'digit')})) then
                        vec ,= p
                fi
        od;
        Vector([vec[]], 'datatype' = prime)
end:
> 

a:=CodeTools:-Usage(A161786__1(10^6,10^6)):numelems(a);

memory used=3.52GiB, alloc change=3.26MiB, cpu time=30.59s, real time=30.63s, gc time=2.52s

42720

> 

CodeTools:-Profiling:-Profile(A161786__1);
numelems(A161786__1(10^6,10^6));

42720

> 

CodeTools:-Profiling:-PrintProfiles(A161786__1);

A161786__1
A161786__1 := proc(m::nonnegint, n::posint := 1, $)
local p, vec;
     |Calls Seconds  Words|
PROC |    1  39.531 475133042|
   1 |    1   0.000     13| p := ifelse(n = 1,1,ithprime(n-1));
   2 |    1   0.000    122| vec := DEQueue();
   3 |    1   1.046      7| to m do
   4 |1000000  37.552 466278938|     if member(4,rhs~({StringTools['
                                  CharacterFrequencies'](nprintf("%d",
                                  (p := nextprime(p))),'digit')})) then
   5 |42720   0.277 1754551|         vec ,= p
                                end if
                            end do;
   6 |    1   0.656 7099411| Vector([vec[]],('datatype') = prime)
end proc
 

> 

 

Download Aseq.mw

avoid large exp...

dharr 8270
October 09 2023
4 1

(worksheet not displaying right now)

contours.mw

Replace G0 with

G0:=proc(y) if not type(y,numeric) then return 'procname'(args) end if;
    if y<0 then exp(y)/(1+exp(y)) else 1/(exp(-y)+1) end if;
   end proc;

vector copy...

dharr 8270
October 08 2023
2 5

Use prev[1..-1] := curr[1..-1]; to copy the contents regardless of size (or the for loop), rather than prev:=curr; You cannot assign to a formal parameter.

This has nothing to do with the size - your first attempt had prev=curr, not prev:=curr so did not work either.

one way...

dharr 8270
October 08 2023
2 11

Here's one way to do it. Of course you could do it all in one loop.

> 

restart;
> 

with(ImageTools):

Import some prn files plt1.prn, plt2.prn,plt3.prn

> 

for i to 3 do
  img[i]:=Read(cat(interface(worksheetdir),"/plt",i,".png"));
end do:

Images are 600x800 pixels

> 

lowerbound(img[1]);
upperbound(img[1]);

1, 1, 1

600, 800, 3

Pad 50 pixels white all around

> 

for i to 3 do
  padded[i]:=PadImage(img[i],border=50,fill=1,reindex=true);
end do:
> 

lowerbound(padded[1]);
upperbound(padded[1]);

1, 1, 1

700, 900, 3

Write out padded files as .tiff

> 

for i to 3 do
  Write(cat(interface(worksheetdir),"/plt",i,".tiff"),padded[i]);
end do;

11270

11284

11284

> 

NULL

Download framepng.mw

Test_prn_files.zip

singularity...

dharr 8270
October 05 2023
1 1

if the singularity is replaced by the limit of the function at the origin, then you can make some progress, though I'm being careless about checking consistency of different limit directions. What do we know about the signs of the parameters?

> 

restart
> 

R0 := -(1-tanh((1/2)*l*sqrt(k*y^2+x^2)*d)^2)*x/(k*y^2+x^2)+2*tanh((1/2)*l*sqrt(k*y^2+x^2)*d)*x/(l*(k*y^2+x^2)^(3/2)*d)

-(1-tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)^2)*x/(k*y^2+x^2)+2*tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)*x/(l*(k*y^2+x^2)^(3/2)*d)

> 

eval(R0, {x = 0, y = 0})

Error, numeric exception: division by zero

However, the function seems smooth, aside from the isolated singularity at the origin

> 

params := {d = 1, k = 1, l = 1}; plot3d(eval(R0, params), x = -.1 .. .1, y = -.1 .. .1)

{d = 1, k = 1, l = 1}

So we should be able to evaluate the limit.

> 

limit(R0, {x = 0, y = 0})

limit(-(1-tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)^2)*x/(k*y^2+x^2)+2*tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)*x/(l*(k*y^2+x^2)^(3/2)*d), {x = 0, y = 0})

Aagh. But

> 

limit(R0, x = 0); `assuming`([limit(R0, y = 0)], [x::real]); limit(%, x = 0)

0

(-4*exp(abs(x)*l*d)*l*d*abs(x)+2*exp(2*abs(x)*l*d)-2)/(x*l*d*abs(x)*(exp(2*abs(x)*l*d)+2*exp(abs(x)*l*d)+1))

0

We could proceed slowly working out all the limits of the required derivatives, e.g.

> 

diff(R0, x); limit(%, x = 0); `assuming`([limit(%, y = 0)], [y::real])

tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)*l*d*x^2*(1-tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)^2)/(k*y^2+x^2)^(3/2)-(1-tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)^2)/(k*y^2+x^2)+3*(1-tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)^2)*x^2/(k*y^2+x^2)^2+2*tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)/(l*(k*y^2+x^2)^(3/2)*d)-6*tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)*x^2/(l*(k*y^2+x^2)^(5/2)*d)

(-4*exp((k*y^2)^(1/2)*l*d)*(k*y^2)^(7/2)*l*d+2*k^3*y^6*exp(2*(k*y^2)^(1/2)*l*d)-2*k^3*y^6)/(k*y^2*(k*y^2)^(7/2)*l*d*(exp((k*y^2)^(1/2)*l*d))^2+2*k*y^2*(k*y^2)^(7/2)*l*d*exp((k*y^2)^(1/2)*l*d)+k*y^2*(k*y^2)^(7/2)*l*d)

(1/6)*d^2*l^2

> 

diff(R0, y); limit(%, x = 0)

tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)*l*d*k*y*(1-tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)^2)*x/(k*y^2+x^2)^(3/2)+3*(1-tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)^2)*x*k*y/(k*y^2+x^2)^2-6*tanh((1/2)*l*(k*y^2+x^2)^(1/2)*d)*x*k*y/(l*(k*y^2+x^2)^(5/2)*d)

0

> 

NULL

Download mtaylor.mw

solve/identity...

dharr 8270
October 02 2023
3 2


 

> 

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

profiling...

dharr 8270
October 01 2023
2 0

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)

Mathematical Expression...

dharr 8270
September 25 2023
3 1

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

Conjugates...

dharr 8270
September 18 2023
0 0

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

unapply...

dharr 8270
September 16 2023
0 0

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

permutation...

dharr 8270
September 16 2023
1 0

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

many solutions...

dharr 8270
September 14 2023
3 3

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]

correct solution...

dharr 8270
September 12 2023
1 0

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

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