dharr

Dr. David Harrington

6120 Reputation

21 Badges

19 years, 330 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 professor of chemistry at the University of Victoria, BC, Canada, where 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

solve, identity will automatically do all the work for you

solve(identity(E1, xi), {k, lambda, w, A[0], A[1], B[1]})

See the attached worksheet for the results.

solve_identity.mw

The pseudo inverse of A seems to work here, but off-hand I can't justify this.

Groebner.mw

 

Data is in the startup code as a hack since you didn't give the datafile. I used fsolve and directly in p__eng1, but NonlinearFit still doesn't return in reasonable time, though you might want to wait longer than I did. As usual in nonlinearfitting, a good initial guess at the solution is probably needed, and reducing the ranges over which you search, if you know them.

test.mw

 

You can remove the solutions you don't want with

remove(has,{COEFFS},[A[0]=0, A[1]=0, B[1]=0]);

If you want to get rid of the RootOfs, add explicit to solve. But if you know the signs of the other parameters add them after the solve, e.g., solve(...) assuming k>0;

eqns.mw

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)

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

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

 

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

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

You may want to add further assumptions.

Download aa.mw

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

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

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

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

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

1 2 3 4 5 6 7 Last Page 1 of 64