## dharr

Dr. David Harrington

## 5116 Reputation

19 years, 101 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.

## integer solutions using algcurves...

Find the smallest positive solution to eq below.

 >
 >
 >

Smallest (!) positive solution is

 >

We more-or-less follow the method given in https://www.agftutoring.com/x-yz-y-xz-z-xy-4/
Rearrange as a polynomial

 >

Convert this homogeneous polynomial to a bivariate polynomial - elliptic curves have genus = 1

 >

The j_invariant has the same value despite the various transformations (birational equivalence)

 >

Maple finds a Weierstrass form

 >

We would like to have it in the standard form , with A and B integers.

Try simply scaling X and Y and the whole expression; we also need to change the sign of X.

 >

 >

Conversion between f and w

 >

Test with a known solution (can be found using method below)

 >

Conversion between w and f

 >

Find the discriminant of the cubic

 >

To find the first solution, we want either , or a y value for which  divides the disciminant, for which the corresponding x is an integer  (Nagell-Lutz) .
So we test all the divisors of

 >
 >

We found two solutions - convert back to a,b,c form and check. First one is OK for eq2, but not for the original eq because a denominator is zero

 >

Error, numeric exception: division by zero

So the only solution is the second one. (This is the single generator of the torsion group.)

 >

To find other points we combine existing point(s) until we find a positive one
Define addition of points in the group algebra - see Wikipedia

A is the coefficient in y^2=x^3+A*x+B

 > algadd := proc(P::[rational,rational], Q::[rational,rational], A::rational)   local s,xP,yP,xQ,yQ,xR,yR;   xP,yP := P[];   xQ,yQ := Q[];   if xP <> xQ then     s := (yP - yQ)/(xP - xQ);     xR := s^2 - xP - xQ;     yR := yP - s*(xP - xR);     [xR, -yR];   elif yP = -yQ then     [0, 0] # infinity point not representable   else # P = Q,sum is tangent     s := (3*xP^2 + A)/(2*yP);     xR := s^2 - 2*xP;     yR := yP - s*(xP - xR);     [xR, -yR]   end if; end proc:
 >

 >

Incrementally add P[1] unitil we get a positive solution

 >
 >

 >

 >

## dcoeffs...

dcoeffs can do this - the order is determined by the type of derivative, not the way the ode is written; they are the same in your example so hope this generalizes the way you want.

 >
 >
 >

 >

 >

## tiny values...

 >
 > with(Student[MultivariateCalculus]):
 >
 >

Values are close to zero, so need to scale function to reasonable values

 >

 >
 >
 >

 >

## LinearSolve...

This can be done with LinearSolve as shown. c1 etc are very complicated expressions that perhaps can be simplified, but I did not wait to see if that were the case. This method assumes that the equations are independent (Determinant(A)<>0), and pays no attention to any relationships that might exist between the Bessel functions.

 > restart;
 > with(LinearAlgebra):
 > fn1 := (2*BesselK(1, alpha1*sigma)*T1*alpha1^3*b1*c+2*BesselK(1, alpha1*sigma)*T1*T2*alpha1*b1*c+BesselK(1, alpha1*sigma)*alpha1*b1*c+BesselK(1, alpha1*sigma)*alpha1*b1+BesselK(0, alpha1*sigma))*c1+(-2*BesselI(1, alpha1*sigma)*T1*alpha1^3*b1*c-2*BesselI(1, alpha1*sigma)*T1*T2*alpha1*b1*c-BesselI(1, alpha1*sigma)*alpha1*b1*c-BesselI(1, alpha1*sigma)*alpha1*b1+BesselI(0, alpha1*sigma))*c2+(2*BesselK(1, alpha2*sigma)*T1*alpha2^3*b1*c+2*BesselK(1, alpha2*sigma)*T1*T2*alpha2*b1*c+BesselK(1, alpha2*sigma)*alpha2*b1*c+BesselK(1, alpha2*sigma)*alpha2*b1+BesselK(0, alpha2*sigma))*c3+(-2*BesselI(1, alpha2*sigma)*T1*alpha2^3*b1*c-2*BesselI(1, alpha2*sigma)*T1*T2*alpha2*b1*c-BesselI(1, alpha2*sigma)*alpha2*b1*c-BesselI(1, alpha2*sigma)*alpha2*b1+BesselI(0, alpha2*sigma))*c4+ur/alpha^2+(-2*A1*K*BesselI(1, k*sigma)*T1*b1*c*k^3+2*B1*K*BesselK(1, k*sigma)*T1*b1*c*k^3-2*A1*K*BesselI(1, k*sigma)*T1*T2*b1*c*k+2*B1*K*BesselK(1, k*sigma)*T1*T2*b1*c*k-2*A1*BesselI(1, k*sigma)*T1*b1*c*k*k1+2*B1*BesselK(1, k*sigma)*T1*b1*c*k*k1-A1*K*BesselI(1, k*sigma)*b1*c*k+B1*K*BesselK(1, k*sigma)*b1*c*k-A1*K*BesselI(1, k*sigma)*b1*k+B1*K*BesselK(1, k*sigma)*b1*k+A1*K*BesselI(0, k*sigma)+B1*K*BesselK(0, k*sigma))*EZ;

 > fn2 := (-2*BesselK(1, alpha1)*T1*alpha1^3*b1*c-2*BesselK(1, alpha1)*T1*T2*alpha1*b1*c-BesselK(1, alpha1)*alpha1*b1*c-BesselK(1, alpha1)*alpha1*b1+BesselK(0, alpha1))*c1+(2*BesselI(1, alpha1)*T1*alpha1^3*b1*c+2*BesselI(1, alpha1)*T1*T2*alpha1*b1*c+BesselI(1, alpha1)*alpha1*b1*c+BesselI(1, alpha1)*alpha1*b1+BesselI(0, alpha1))*c2+(-2*BesselK(1, alpha2)*T1*alpha2^3*b1*c-2*BesselK(1, alpha2)*T1*T2*alpha2*b1*c-BesselK(1, alpha2)*alpha2*b1*c-BesselK(1, alpha2)*alpha2*b1+BesselK(0, alpha2))*c3+(2*BesselI(1, alpha2)*T1*alpha2^3*b1*c+2*BesselI(1, alpha2)*T1*T2*alpha2*b1*c+BesselI(1, alpha2)*alpha2*b1*c+BesselI(1, alpha2)*alpha2*b1+BesselI(0, alpha2))*c4+ur/alpha^2+(2*A1*K*BesselI(1, k)*T1*b1*c*k^3-2*B1*K*BesselK(1, k)*T1*b1*c*k^3+2*A1*K*BesselI(1, k)*T1*T2*b1*c*k-2*B1*K*BesselK(1, k)*T1*T2*b1*c*k+2*A1*BesselI(1, k)*T1*b1*c*k*k1-2*B1*BesselK(1, k)*T1*b1*c*k*k1+A1*K*BesselI(1, k)*b1*c*k-B1*K*BesselK(1, k)*b1*c*k+A1*K*BesselI(1, k)*b1*k-B1*K*BesselK(1, k)*b1*k+A1*K*BesselI(0, k)+B1*K*BesselK(0, k))*EZ;

 > fn3 := -T1*alpha1*(alpha1^2+T2)*(BesselK(0, alpha1*sigma)*alpha1*b2*sigma+BesselK(1, alpha1*sigma)*b2*s1+BesselK(1, alpha1*sigma)*b2+BesselK(1, alpha1*sigma)*sigma)*c1/sigma+T1*alpha1*(alpha1^2+T2)*(-BesselI(0, alpha1*sigma)*alpha1*b2*sigma+BesselI(1, alpha1*sigma)*b2*s1+BesselI(1, alpha1*sigma)*b2+BesselI(1, alpha1*sigma)*sigma)*c2/sigma-T1*alpha2*(alpha2^2+T2)*(BesselK(0, alpha2*sigma)*alpha2*b2*sigma+BesselK(1, alpha2*sigma)*b2*s1+BesselK(1, alpha2*sigma)*b2+BesselK(1, alpha2*sigma)*sigma)*c3/sigma+T1*alpha2*(alpha2^2+T2)*(-BesselI(0, alpha2*sigma)*alpha2*b2*sigma+BesselI(1, alpha2*sigma)*b2*s1+BesselI(1, alpha2*sigma)*b2+BesselI(1, alpha2*sigma)*sigma)*c4/sigma+T1*EZ*k*(K*k^2+K*T2+k1)*(-A1*BesselI(0, k*sigma)*b2*k*sigma-B1*BesselK(0, k*sigma)*b2*k*sigma+A1*BesselI(1, k*sigma)*b2*s1-B1*BesselK(1, k*sigma)*b2*s1+A1*BesselI(1, k*sigma)*b2+A1*BesselI(1, k*sigma)*sigma-B1*BesselK(1, k*sigma)*b2-B1*BesselK(1, k*sigma)*sigma)/sigma;

 > fn4 := T1*alpha1*(alpha1^2+T2)*(BesselK(1, alpha1)*b2*s1+BesselK(0, alpha1)*alpha1*b2+BesselK(1, alpha1)*b2-BesselK(1, alpha1))*c1-T1*alpha1*(alpha1^2+T2)*(BesselI(1, alpha1)*b2*s1-BesselI(0, alpha1)*alpha1*b2+BesselI(1, alpha1)*b2-BesselI(1, alpha1))*c2+T1*alpha2*(alpha2^2+T2)*(BesselK(1, alpha2)*b2*s1+BesselK(0, alpha2)*alpha2*b2+BesselK(1, alpha2)*b2-BesselK(1, alpha2))*c3-T1*alpha2*(alpha2^2+T2)*(BesselI(1, alpha2)*b2*s1-BesselI(0, alpha2)*alpha2*b2+BesselI(1, alpha2)*b2-BesselI(1, alpha2))*c4-T1*EZ*k*(K*k^2+K*T2+k1)*(A1*BesselI(1, k)*b2*s1-A1*BesselI(0, k)*b2*k-B1*BesselK(1, k)*b2*s1-B1*BesselK(0, k)*b2*k+A1*BesselI(1, k)*b2-B1*BesselK(1, k)*b2-A1*BesselI(1, k)+B1*BesselK(1, k));

Eqns are linear in c1,c2,c3,c4 - convert to matrix form

 > A,b:=GenerateMatrix([fn1,fn2,fn3,fn4],[c1,c2,c3,c4]):

Solve the system with simpler variables

 > AA:=Matrix(4,4,symbol=x):bb:=Vector(4,symbol=y):
 > CC:=LinearSolve(AA,bb):

Substitute the more complicated expressions into the solution

 > C:=eval(CC,{entries(AA=~A,'nolist'),entries(bb=~b,'nolist')}):
 > soln:={c1=C[1],c2=C[2],c3=C[3],c4=C[4]}:

Check the solution

 > simplify(eval([fn1,fn2,fn3,fn4],soln));

## live connection...

I assume your tools->options->network has "enable maple cloud connection" checked? (becomes active next time you open Maple)

## intentional...

I agree that the documentation suggests that Q and rp are the same size "r". However, a look at the code using showstat(RowEchelonTransform) suggests that the length of Q (=P in the code) need not be the same as rp. P is initialized as a 1..n row vector (n=number of rows presumably) (line 7), but on line 9 only the entries from 1 to min(r,n-1) are output. Presumably r is the rank and in your case min(16,15) = 15.

Since the length of Q is the number of swaps needed to construct P, perhaps you don't need to do the last one if it is full rank, just do all the swaps in Q. I suspect this will work, but if not you could just do RowReduce as a workaround if you only want the row reduced form.

Edit: I confirmed this works for the reduced form (redflag = true). The redflag = false case gives different answers for RowEchelonTransform and RowReduce, but the nonreduced form is not unique. So I think the swap procedure is correct.

I submitted an SCR for the help page.

 > restart;
 > with(LinearAlgebra[Modular]):interface(rtablesize=39):
 > p:=29:
 > M1:=[[-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,1,-1,-1,-1,-1],[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,-1,-1,1,1,1,1],[0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,-1,0,0,0,1],[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,-1],[0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0],[0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,-1,0,0,1,0,0,0,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0]]:
 > nrows:=nops(M1);ncols:=nops(M1[1]);

 > M:=Mod(p,M1,integer[]):

RowEchelonTransform - U1 is altered in place.

 > U1:=Copy(p,M): Q,rp,d:=RowEchelonTransform(p, U1, true, true, true, false):
 > Q;nQ:=numelems(Q); rp;r:=numelems(rp);

 > Rank(p,M); RankProfile(p,M);

Construct P by doing the row swaps specified in Q

 > P:=Create(p,nrows,nrows,identity,integer[]): for i to nQ do   Swap(p,P,i,P,Q[i]); end do:

Construct U from U1 - first r columns from U1, rest are identity (in this case U is just the identity)

 > U:=Create(p,nrows,nrows,identity,integer[]): U[..,1..r]:=U1[..,1..r]:

Row reduced form is U.P.M

 > RR:=Multiply(p,U,Multiply(p,P,M)):

Do the same thing with RowReduce

 > RR2:=Copy(p,M): RowReduce(p,RR2,nrows,ncols,ncols,'det',0,'rank',0,0,true);
 > RR2,det,rank:
 > EqualEntries(RR,RR2);

 >

## difficult to track down...

Your side note question is the easiest - the spaces are missing commas. The code is something like this

```P:=proc();
userinfo(1,P,`entered with args:`,args);
7;
end proc:
infolevel[P]:=1;
P(as,3,5);```

which gives

P: entered with args: as 3 5

For your main question, I don't think there is a simple answer - many internal routines are called and internal tables probably retain some information. For your first example, the message

Solving each unknown as a function of the next ones using the order: [y(t), x(t)]

comes from `PDEtools/sdsolve`, but I don't see where the next message comes from. You could set printlevel to a high value to track it down, but even if you did, the solution (perhaps forget on a deep routine?) would have to be different for different odes.

(I encountered the same problem with SolveTools:-Polynomials system and infolevel[solve].)

## output=basis...

output=basis is for exact solutions; earlier you had integer parameters and an analytical solution in terms of Heun functions was possible (not sure if this is the same equation). So you could try leaving those parameters as variables, find the exact solution, and then plug the values of the parameters into the solutions. I gave up on this after some time computing, but it might ultimately work. (The general solution is then a*(solution 1)+b*(solution 2), where a and b depend on the initial conditions.) Otherwise, for a numerical solution you need to choose initial conditions, as you succeeded to do in an earlier case.

## product only...

My interpretation is that convert(...,unit_free) only removes the units from a product of a number and unit. So it removed the 1/um^2 but not the unit in the argument of the exponential ("data structure with units embedded in it").

`subsindets(expr, anything, convert, unit_free)`

will remove all the units.

## missing function argument...

You had diff(E,psi), which evaluates to 0 instead of diff(E(psi),psi). Note also gamma has a special meaning unless you use "local gamma". Changes in red. I didn't proceed beyond the mess you now get for fin1.

Note: Maple's JacobiSN is spelled with J not j, and has two arguments, so you will have to fix that.

0123.mw

## Matrix-based solution...

@MaPal93 sent me some notes about how argmin was derived, so I could proceed with a matrix-based method, similar to in the earlier thread. This can probably be refined for the second part of the problem, but I used it in the maximization part to avoid ever creating a problematic denominator, thus avoiding the issue. For the calibration problem, there are a few restrictions on lambda__2, then solve a quadratic for lambda__3, then a linear equation for lambda__1. So finding conditions for positive solutions should not be intractable, though perhaps tedious.

MatrixNew2.mw

Edit: There are no positive solutions for the calibration problem (meaning the conjecture is not correct), and a more detailed analysis is given below.

## probably not possible...

Unless there is some dramatic cancelling of terms, then probably not. A 32x32 determinant has 32! terms. Now 32!/16! = 12576278705767096320000 = 1.26 x 10^22, which is how many times longer the calculation for a 32x32 determinant will be over a 16x16 one as an estimate.

On the other hand, your matrix has a lot of structure, which you presumably know about in detail. Did the 16x16 one have 16! terms or was it much simpler? Do you have a guess as to what the answer is, or what form it takes? Perhaps it can be divided into simpler tasks (Determinants of submatrices, perhaps). Do you really want this determinant, or is it just a step to something else that can be calculated without it (solutions to a linear system for example)? Do you just need some information about it? - zero or non-zero perhaps? Do you know anything about the signs of the variables?

Without some guidance about the nature of the problem, which might suggest a clever solution, I'd say you are out of luck.

## solve vs SolveTools:-PolynomialSystem...

This partly duplicates @mmcdara's analysis (I voted up!) that used numerical values for the parameters, but with the analytical solutions, and shows how to interpret the PolynomialSystem answer and do the backsubstitution (part of your second question). In this case solve is simpler (and presumably the same as if backsubstitute=true was chosen, but I didn't check), but had there been restrictions on lambda__2 and lambda__3, then fishing out the required solutions from all of them probably would have been easier starting from the PolynomialSystem, backsubstitute=false form, as was the case in your earlier problem.

 > restart;

@mmcdara's solve solution (without assumptions on the variables).

The interpretation is, as @mmcdara noted, that one can choose any lambda__2 and lambda__3 and from them find an acceptable lambda__1

 > sol_solve := {lambda__1 = -(sigma__v^2*(p - 1) - sigma__e^2)*(p*gamma*(lambda__2 + lambda__3)*sigma__v^4 + (sigma__e^2*(lambda__2 + lambda__3)*gamma + 2*lambda__2*lambda__3)*sigma__v^2 + 2*lambda__2*lambda__3*sigma__e^2)*sigma__v^2*gamma/(p*gamma^2*(p - 1)*sigma__v^8 - 2*(1/2*gamma*sigma__e^2 + p*lambda__3 + lambda__2)*gamma*sigma__v^6 + (-gamma^2*sigma__e^4 - 2*sigma__e^2*((p + 1)*lambda__3 + 3*lambda__2)*gamma - 4*lambda__2*lambda__3)*sigma__v^4 - 4*sigma__e^2*(sigma__e^2*(lambda__2 + 1/2*lambda__3)*gamma + 2*lambda__2*lambda__3)*sigma__v^2 - 4*lambda__2*lambda__3*sigma__e^4), lambda__2 = lambda__2, lambda__3 = lambda__3}

Extract the numerator and denominator pieces of the lambda__1 solution. It is implicit in this solution that the lambda__2 and lambda__3 chosen must not make the denominator zero.

 > solnm1:=eval(lambda__1,sol_solve): n1=numer(%);d1:=denom(%%);

@MaPal93's solution from SolveTools:-PolynomialSystem using the variable order lambda__1 > lambda__3 > lambda__2 (see the userinfo output)

 > Solution := [[(((-4*sigma__e^4 - 8*sigma__e^2*sigma__v^2 - 4*sigma__v^4)*lambda__2 - 2*gamma*p*sigma__e^2*sigma__v^4 - 2*gamma*p*sigma__v^6 - 2*gamma*sigma__e^4*sigma__v^2 - 2*gamma*sigma__e^2*sigma__v^4)*lambda__3 + (-4*gamma*sigma__e^4*sigma__v^2 - 6*gamma*sigma__e^2*sigma__v^4 - 2*gamma*sigma__v^6)*lambda__2 + gamma^2*p^2*sigma__v^8 - gamma^2*p*sigma__v^8 - gamma^2*sigma__e^4*sigma__v^4 - gamma^2*sigma__e^2*sigma__v^6)*lambda__1 + ((2*gamma*p*sigma__e^2*sigma__v^4 + 2*gamma*p*sigma__v^6 - 2*gamma*sigma__e^4*sigma__v^2 - 4*gamma*sigma__e^2*sigma__v^4 - 2*gamma*sigma__v^6)*lambda__2 + gamma^2*p^2*sigma__v^8 - gamma^2*p*sigma__v^8 - gamma^2*sigma__e^4*sigma__v^4 - gamma^2*sigma__e^2*sigma__v^6)*lambda__3 + (gamma^2*p^2*sigma__v^8 - gamma^2*p*sigma__v^8 - gamma^2*sigma__e^4*sigma__v^4 - gamma^2*sigma__e^2*sigma__v^6)*lambda__2], {((-4*sigma__e^4 - 8*sigma__e^2*sigma__v^2 - 4*sigma__v^4)*lambda__2 - 2*gamma*p*sigma__e^2*sigma__v^4 - 2*gamma*p*sigma__v^6 - 2*gamma*sigma__e^4*sigma__v^2 - 2*gamma*sigma__e^2*sigma__v^4)*lambda__3 + (-4*gamma*sigma__e^4*sigma__v^2 - 6*gamma*sigma__e^2*sigma__v^4 - 2*gamma*sigma__v^6)*lambda__2 + gamma^2*p^2*sigma__v^8 - gamma^2*p*sigma__v^8 - gamma^2*sigma__e^4*sigma__v^4 - gamma^2*sigma__e^2*sigma__v^6 <> 0}]

Extract the solution and the left-hand side of the condition<>0.

 > soln:=Solution[][1][]; cond:=lhs(Solution[][2][]);

The output from the solver with backsubstitute=false is usually three (partial) solutions, which using the solver order here would contain (lambda__1, lambda__2, lambda__3),  (lambda__2, lambda__3) and (lambda__2). The backsubstitution process consists of finding the solutions to the third equation for lambda__2 substituting into the second equation to find the values of lambda__3 and then substituting the values for lambda__3 and lambda__2 into the first solution to find the values of lambda__1. In this case, the lack of the second and third solutions is just a way of saying that lambda__2 and lambda__3 can have any values. Then we solve for lambda__1 as follows

 > solnm2:=solve(soln,lambda__1);

This is the same as the solution solve gave, and the condition is just that the denominator is not zero

 > simplify(solnm1-solnm2); simplify(cond-d1);

 >