Advanced Search

dharr

Dr. David Harrington

8482 Reputation

22 Badges

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

double integral - alternate form...

dharr 8482
July 11 2024
3 0 
int(F[1](x,y),[x=0..1.,y=0..1.]);

returns 0.2080471649 (This is Maple 2024, which also has a problem with the nested form)

couple of things...

dharr 8482
July 11 2024
1 0

You eq (5) label refers to length(%), not the collected expression, so delete length(%). Change coeff(.., lambda^0) to coeff(..,lambda,0).

But your collect is on a ratio. If you want to just work with the numerator, that can be done, as on the attached.

coeff.mw

plots:-inequal...

dharr 8482
July 10 2024
1 1

The following works, though it is slow.

plots:-inequal({-1 < lambda1, -1 < lambda2, lambda1 < 1, lambda2 < 1}, v = -5 .. 5, z = -5 .. 5)


The colored regions are where both eigenvalues are between -1 and 1. Depending on the ranges you want for v and z, you can refine this, and there are other options for inequal that might increase resolution or efficiency.

2d_implicit_plot_[v_z].mw

 

hyperlink...

dharr 8482
July 07 2024
1 1

For FriendshipGraph in SpecialGraphs in GraphTheory GraphTheory[SpecialGraphs][FriendshipGraph] works.

RootOf reliability...

dharr 8482
July 06 2024
7 5

I'm not sure if this is an answer to the specific question, but here is my take on it. Basically, RootOf is very reliable for roots of polynomials, but if there are square roots or cube roots inside it, it is better IMO to reformulate it without these. Here is a way to do that for this case.

dsolve without method

> 

ode:=diff(y(x), x) = (3*x - y(x) + 1)/(3*y(x) - x + 5);
ic:=y(0)=0;
dsolve({ode,ic},implicit);
isol:=eval(lhs(%),y(x)=y);

diff(y(x), x) = (3*x-y(x)+1)/(3*y(x)-x+5)

y(0) = 0

-(2/3)*ln(-(3+y(x)+x)/(x+1))-(1/3)*ln((-1-y(x)+x)/(x+1))-ln(x+1)+(2/3)*ln(3)+I*Pi = 0

-(2/3)*ln(-(3+y+x)/(x+1))-(1/3)*ln((-1-y+x)/(x+1))-ln(x+1)+(2/3)*ln(3)+I*Pi

We need to combine the two arguments of the ln to solve this in the first 2 terms, but the 1/3 powers are a problem

> 

combine(isol,ln);
isol2:=combine(-3*isol,ln);

ln(1/(-(3+y+x)/(x+1))^(2/3))+ln(1/((-1-y+x)/(x+1))^(1/3))-ln((1/3)*(x+1)*3^(1/3))+I*Pi

2*ln(-(3+y+x)/(x+1))+3*ln(x+1)+ln((1/9)*(-1-y+x)/(x+1))-(3*I)*Pi

> 

with_y,no_y:=selectremove(has,isol2,y);

2*ln(-(3+y+x)/(x+1))+ln((1/9)*(-1-y+x)/(x+1)), 3*ln(x+1)-(3*I)*Pi

> 

p1:=simplify(exp(with_y-ln(a)));
p2:=simplify(exp(no_y+ln(a)));

(1/9)*(3+y+x)^2*(-1-y+x)/((x+1)^3*a)

-(x+1)^3*a

So now we have p1*p2=exp(0)

> 

eq:=9*(1-p1*p2);

(-1-y+x)*(3+y+x)^2+9

> 

plots:-implicitplot(eq,x=-10..10,y=-2..9);

solve(.., parametric) suggests that except for the case of x=-1, we have the three solutions of the cubic, the same as we would get just with solve.

> 

solve(eq,y,parametric);

piecewise(x+1 = 0, %SolveTools[Engine]({-y^3-6*y^2-12*y+1}, {y}), x+1 <> 0, [{y = (1/6)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(6*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(1/3)*x-7/3}, {y = -(1/12)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(3*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(1/3)*x-7/3+I*sqrt(3)*((1/6)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(6*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3))*(1/2)}, {y = -(1/12)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(3*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(1/3)*x-7/3-I*sqrt(3)*((1/6)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(6*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3))*(1/2)}])

But this seems to miss the real solution for x <~-3 from the plot above - in radical form there are still issues with which square root or which cube root to take

> 

plot([solve(eq,y)],x=-10..10,color=[red, blue, black]);

indexed RootOf is reliable for polynomials

> 

plot([seq(RootOf(eq,y,index=i),i=1..3)],x=-10..10,color=[red, blue, black]);

> 

NULL

 

NULL

Download RootOf_stuff.mw

assuming...

dharr 8482
June 30 2024
1 0

You may want to add further assumptions.

Download aa.mw

Straightforward...

dharr 8482
June 25 2024
2 3

Now that you have corrected the integral with the exp (though one was ecp), the integral is evaluated without any special treatment.

(I had to guess at a value for x.)

> 

restart;
> 

local D; # avoid clash with differential operator

D

> 

Wminus := M*sqrt(Pi)/(2*A*a)*(x*erf(a*x) + (1/(a*sqrt(Pi)))*exp(-a^2*x^2)) + (M/(A*a*sqrt(Pi)))*(Int(exp(-epsilon^2/(4*a^2) - tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (2*m/Pi)*(Int(exp(-tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (M/(A*a*sqrt(Pi)))*(Int((1/(2*a^2) + 3*tau*epsilon)*exp(-epsilon^2/(4*a^2) - tau*epsilon^3), epsilon = 0..infinity)) - (2*m/Pi)*tau^(1/3)*GAMMA(2/3) -m*x;

(1/2)*M*Pi^(1/2)*(x*erf(a*x)+exp(-a^2*x^2)/(a*Pi^(1/2)))/(A*a)+M*(Int(exp(-(1/4)*epsilon^2/a^2-tau*epsilon^3)*(cos(epsilon*x)-1)/epsilon^2, epsilon = 0 .. infinity))/(A*a*Pi^(1/2))+2*m*(Int(exp(-tau*epsilon^3)*(cos(epsilon*x)-1)/epsilon^2, epsilon = 0 .. infinity))/Pi+M*(Int(((1/2)/a^2+3*tau*epsilon)*exp(-(1/4)*epsilon^2/a^2-tau*epsilon^3), epsilon = 0 .. infinity))/(A*a*Pi^(1/2))-2*m*tau^(1/3)*GAMMA(2/3)/Pi-m*x

> 

Wplus := M*sqrt(Pi)/(2*A*a)*(x*erf(a*x) + (1/(a*sqrt(Pi)))*exp(-a^2*x^2)) + (M/(A*a*sqrt(Pi)))*evalf(Int(exp(-epsilon^2/(4*a^2) - tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (2*m/Pi)*evalf(Int(exp(-tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (M/(A*a*sqrt(Pi)))*evalf(Int((1/(2*a^2) + 3*tau*epsilon)*exp(-epsilon^2/(4*a^2) - tau*epsilon^3), epsilon = 0..infinity)) - (2*m/Pi)*tau^(1/3)*GAMMA(2/3) + m*x;

(1/2)*M*Pi^(1/2)*(x*erf(a*x)+exp(-a^2*x^2)/(a*Pi^(1/2)))/(A*a)+M*(Int(exp(-.2500000000*epsilon^2/a^2-1.*tau*epsilon^3)*(cos(epsilon*x)-1.)/epsilon^2, epsilon = 0. .. Float(infinity)))/(A*a*Pi^(1/2))+2*m*(Int(exp(-1.*tau*epsilon^3)*(cos(epsilon*x)-1.)/epsilon^2, epsilon = 0. .. Float(infinity)))/Pi+M*(Int((.5000000000/a^2+3.*tau*epsilon)*exp(-.2500000000*epsilon^2/a^2-1.*tau*epsilon^3), epsilon = 0. .. Float(infinity)))/(A*a*Pi^(1/2))-2*m*tau^(1/3)*GAMMA(2/3)/Pi+m*x

> 

params:= {x=1e-3, M = 6.3*10^17, t = 7200, a = 10^3, m = 1, D = 1.0*10^(-5), Omega = 1.7*10^23, tau = 1.0*10^(-14)}:
params:= {A = eval(tau/(D*Omega*t),params)} union params;
> 

 

{A = 0.8169934640e-36, D = 0.1000000000e-4, M = 0.6300000000e18, Omega = 0.1700000000e24, a = 1000, m = 1, t = 7200, tau = 0.1000000000e-13, x = 0.1e-2}

> 

evalf(eval(Wminus,params));

0.7711380641e48

> 

evalf(eval(Wplus,params));

0.7711380641e48

> 

NULL

Download integration.mw

solution...

dharr 8482
June 23 2024
3 1

I agree with @nm; isolve is frequently weak, especially here where the equations are nonlinear. Here is another ad-hoc way that works.

> 

restart;
> 

eq1 := a*b + c = 2020;
eq2 := b*c + a = 2021;
 

a*b+c = 2020

b*c+a = 2021

Eliminate b, then isolve can handle the remaining equation in a and c

> 

eliminate({eq1, eq2}, {b});
beq,eq3:=%[];

[{b = -(c-2020)/a}, {a^2-c^2-2021*a+2020*c}]

{b = -(c-2020)/a}, {a^2-c^2-2021*a+2020*c}

> 

[isolve(eq3)];
eq3solns:=remove(has,%,a=0); #avoid division by zero

[{a = 0, c = 0}, {a = 0, c = 2020}, {a = 673, c = 674}, {a = 673, c = 1346}, {a = 896, c = 900}, {a = 896, c = 1120}, {a = 1125, c = 900}, {a = 1125, c = 1120}, {a = 1348, c = 674}, {a = 1348, c = 1346}, {a = 2021, c = 0}, {a = 2021, c = 2020}]

[{a = 673, c = 674}, {a = 673, c = 1346}, {a = 896, c = 900}, {a = 896, c = 1120}, {a = 1125, c = 900}, {a = 1125, c = 1120}, {a = 1348, c = 674}, {a = 1348, c = 1346}, {a = 2021, c = 0}, {a = 2021, c = 2020}]

Substitute back into the equation for b and select the solutions where b is an integer

> 

map(x->eval(beq,x) union x,eq3solns);
select(x->eval(b,x)::integer,%);

[{a = 673, b = 2, c = 674}, {a = 673, b = 674/673, c = 1346}, {a = 896, b = 5/4, c = 900}, {a = 896, b = 225/224, c = 1120}, {a = 1125, b = 224/225, c = 900}, {a = 1125, b = 4/5, c = 1120}, {a = 1348, b = 673/674, c = 674}, {a = 1348, b = 1/2, c = 1346}, {a = 2021, b = 2020/2021, c = 0}, {a = 2021, b = 0, c = 2020}]

[{a = 673, b = 2, c = 674}, {a = 2021, b = 0, c = 2020}]

> 

NULL

Download isolve.mw

local variable...

dharr 8482
June 22 2024
3 1

If you change local A:=0,B to global A:=0,B it works as expected. When you print out A, you are printing out the local A, and then you use that value to determine the value to be saved.

The help for read says "If the file is in Maple language format, the statements in the file are read and executed as if they were being entered into Maple interactively", which I take to mean that A := 6 assigns a value to global A.

You could adjust by being more explicit about the locals and globals:
 

> 

restart;
> 

currentdir();
try  #remove TMP.m first time
    FileTools:-Remove("TMP.m");
catch:
    NULL;
end try;

#this function is called many times.
foo:=proc()
 local A:=0,B;
 A:=A+10;
 if FileTools:-Exists("TMP.m") then
     print("TMP.m exists, will read its content");
     B:=A; #save copy
     read "TMP.m";
     print("Read old value of A from file",:-A);
     A:= :-A + B; # update running total
     print("Now A is ",A," will now save this new value");     
     save A,"TMP.m";
  else
     print("TMP.m do not exist, first time saving to it");
     save A,"TMP.m";
  fi;
  print("A=",A, :-A);
end proc:
 

"C:\Users\dharr\OneDrive - University of Victoria\Desktop"

> 

foo();

"TMP.m do not exist, first time saving to it"

"A=", 10, A

> 

foo();

"TMP.m exists, will read its content"

"Read old value of A from file", 10

"Now A is ", 20, " will now save this new value"

"A=", 20, 10

> 

foo();

"TMP.m exists, will read its content"

"Read old value of A from file", 20

"Now A is ", 30, " will now save this new value"

"A=", 30, 20

> 

foo();

"TMP.m exists, will read its content"

"Read old value of A from file", 30

"Now A is ", 40, " will now save this new value"

"A=", 40, 30

> 

 


 

Download T.mw

if statement...

dharr 8482
June 19 2024
0 3

Note that sigma(a)*sigma(b) is not the same as sigma(a*b), which it seemed to me you were assuming. I think this answers your question about the if statement.

> 

restart
> 

with(numtheory)
> 

sigma(12)*sigma(6)

336

> 

sigma(12*6)

195

> 

p := sigma(12)-2*12

4

> 

if p>0 then
  print(abundant)
elif p=0 then
  print(perfect)
else
  print(deficient)
end if;

abundant

> 

``

Download abundant_edge.mw

variation...

dharr 8482
June 15 2024
2 4

Here's a variation on @Ronan's answer.

> 

restart;
> 

for n to 4 do
   A[n]:=Matrix(n,n,(i,j)->local m;
                           ifelse(i+j<n+2,
                                  add(f[i-m+j-1]*binomial(i-1,m)*(-1)^m,
                                      m=0..i-1),
                                  0));
end do;

Vector(1, {(1) = f[1]})

Matrix(2, 2, {(1, 1) = f[1], (1, 2) = f[2], (2, 1) = f[2]-f[1], (2, 2) = 0})

Matrix(3, 3, {(1, 1) = f[1], (1, 2) = f[2], (1, 3) = f[3], (2, 1) = f[2]-f[1], (2, 2) = f[3]-f[2], (2, 3) = 0, (3, 1) = f[3]-2*f[2]+f[1], (3, 2) = 0, (3, 3) = 0})

Matrix(%id = 36893628080372195684)

Download binomial.mw

bug reporting...

dharr 8482
June 15 2024
1 0

I agree with your analysis, and that it is a bug. To report it, you can (logged into Mapleprimes) choose the "more" tab, and then "submit software change request". I usually copy the mapleprimes link and paste that in the box where you are asked for any further information.

comments...

dharr 8482
June 04 2024
2 3

I don't know the specific answer to question 1. If you extend to Gamma negative you don't see a region. Anyway, as you say, it gives the correct number of solutions in the different regions shown. Other explanations below.

> 

restart;

Here @acer's analysis

> 

Eq := -8*(rho + 1)^4*Lambda^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*Lambda^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*Lambda^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*Lambda + Gamma^2*(rho + 1)*(rho - 1)^2;

-8*(rho+1)^4*Lambda^4+12*(rho+1)^3*Gamma*(rho-1)*Lambda^3-5*(rho+1)^2*(-4/5+Gamma^2*rho^2+2*(-2/5-Gamma^2)*rho+Gamma^2)*Lambda^2-4*(rho+1)*Gamma*(rho^2-1)*Lambda+Gamma^2*(rho+1)*(rho-1)^2

> 

with(RootFinding:-Parametric):
> 

m := CellDecomposition([Eq=0, Lambda>0, Gamma>0, rho>-1, rho<1], [Lambda]):
> 

CellPlot(m, samplepoints, view=[0..10,-2..2], size=[400,300]);

Number of solutions in regions 1..5

> 

NumberOfSolutions(m)[1..5];

[[1, 0], [2, 1], [3, 1], [4, 0], [5, 0]]

The chosen sample points within the region (as shown on the plot)

> 

m:-SamplePoints;

[[Gamma = 1, rho = -2], [Gamma = 1, rho = 0], [Gamma = 2, rho = 0], [Gamma = 1, rho = 2], [Gamma = 2, rho = 2]]

Value of Lambda at sample point 2

> 

SampleSolutions(m,2);

[[Lambda = .5622851345]]

> 

SampleSolutions(m,3);

[[Lambda = .5161075747]]

The plot3d that I think you are worried about. For the region Lambda>0, there is everywhere a solution (the top sheet), in agreement with the analysis above that shows 1 real solution in regions 2and 3

> 

plot3d([allvalues(RootOf(Eq,Lambda))],rho=-1..1,Gamma=0..5, grid = [20,20]);

> 

NULL

Download cells.mw

 

No...

dharr 8482
May 27 2024
1 5

If you mean multiply by some function of the variables then the answer is no. The answer is somehow "closer" if you change the red k/2 to -k/2. What are the assumption on k? If one is allowed to multiply by matrices, then the answer might be different. (Matrices not showing correctly on Mapleprimes.

> 

restart
> 

with(LinearAlgebra)
> 

L := Matrix([[(1+I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp(-(1/2)*k), I*rho*(exp((1/2)*k)-exp(-(1/2)*k))/lambda], [I*rho*(exp((-k)*(1/2))-exp((1/2)*k))/lambda, (1-I*w*sqrt((1-cos(k))/(1+cos(k)))/lambda)*exp((1/2)*k)]])

Matrix(%id = 36893490028717500884)

> 

X := Matrix([[-I*lambda*(1/2)-(1/2)*w, -rho], [rho, I*lambda*(1/2)+(1/2)*w]])

Matrix(%id = 36893490028738504332)

To get the off diagonals of L to be those of X, we need to multiply by

> 

F := L*lambda/(I*(exp(-(1/2)*k)-exp((1/2)*k)))

Matrix(%id = 36893490028738496500)

> 

F__2 := `assuming`([simplify(map(convert, F, exp))], [k > 0, k < (1/2)*Pi])

Matrix(%id = 36893490028738516380)

Diagonals have different sums, so cannot be the same.

> 

simplify(Trace(F__2)); Trace(X)

(w*(exp(k)-1)*tan((1/2)*k)+I*(exp(k)+1)*lambda)/(exp(k)-1)

0

> 

NULL

Download XintoL2.mw

PolynomialSystem...

dharr 8482
May 26 2024
0 4

SolveTools:-PolynomialSystem with engine=triade can solve this. PolynomialSystems allows choice of different engines, one of which is groebner, so you can experiment if it gets stuck.

solve_test.mw

First 18 19 20 21 22 23 24 Last Page 20 of 83