## dharr

Dr. David Harrington

## 6120 Reputation

19 years, 330 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

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.

## solve,identity...

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

## pseudo-inverse...

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

Groebner.mw

## progress...

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

## remove...

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

## double integral - alternate form...

`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...

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...

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.

## RootOf reliability...

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);

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);

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

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

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

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

 > 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);

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]);

 >

## assuming...

You may want to add further assumptions.

## Straightforward...

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

 > 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;

 > 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;

 > 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;
 >

 > evalf(eval(Wminus,params));

 > evalf(eval(Wplus,params));

 >

## solution...

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;

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

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

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

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,%);

 >

## local variable...

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:

 > foo();

 > foo();

 > foo();

 > foo();

 >

## if statement...

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.

 >
 >
 >

 >

 >

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

 >