dharr

Dr. David Harrington

8235 Reputation

22 Badges

20 years, 341 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 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 replies submitted by dharr

@MaPal93 You had not done the nondimensionalization with dims at the beginning. If you do that the scaling is easy to see. I added code to reduce multiple numerical evaluations since the expression is getting complicated.

deriv_beta1.mw

I didn't follow exactly what you want to do for auto scaling, but since the powers are integers, maybe degree does what you want, applied to the numerators or denominators after normal(.., expanded).

@MaPal93 Here is a way to get nice plots without artifacts. I tried it only for d(lambda__2)/d(sigma__d), where it works well. Yours was scaled by an extra factor of Gamma__1, so they aren't directly comparable. That could be easily changed; the scaling is somwhat arbitrary. It could be made into a procedure, but figuring out the scaling to make it dimensionless was done in an ad-hoc fashion.

In the plots you showed, there always was an OK scaling so the scaled derivative was a function only of Gamma__1 and Gamma__2. But for your derivatives of alpha etc it seems much harder to achieve that and it might not be possible? But if you have achieved that or it is achievable, then this general method should work on those cases.

Efficient evaluation of derivatives of the original (dimensioned) variables.
Strategy: 1. Plot scaled dimensionless versions of these vs Gamma__1 and Gamma__2. Keep all calculations in terms of dimensionles variables as much as possible.

Create a procedure that takes Gamma__1 and Gamma__2 and produces the scaled derivative, which
2. Evaluates Lambda__1(Gamma__1, Gamma__2) numerically only once, extracting the most positive root
3. Use this to evaluate Lambda__2(Gamma__1, Gamma__2)  numerically only once.
4. Evaluate the deriatives D[1](Lambda__1) etc using fdiff
5. Combine 2, 3, and 4 to give required derivative.
Equations to be solved, their nondimensioalization and solution are in the startup code.

restart;

local gamma:

Non-dim and dim variables. Since there are more dim vars than non-dim vars, we allow sigma__d and sigma__v2  (and sigma__d3 if using the 3rd eqn) -see dims. Can scale by these after the differentiation. (Definitions in startup code.)

'nondims'=nondims;
'dims'=dims;

nondims = {Gamma__1 = sigma__d*sigma__v1*gamma, Gamma__2 = sigma__d*sigma__v2*gamma, Lambda__1 = sigma__d*lambda__1/sigma__v1, Lambda__2 = sigma__d*lambda__2/sigma__v1, Lambda__3 = sigma__d3*lambda__3/sigma__v1}

dims = {gamma = Gamma__2/(sigma__d*sigma__v2), lambda__1 = Lambda__1*Gamma__1*sigma__v2/(Gamma__2*sigma__d), lambda__2 = Lambda__2*Gamma__1*sigma__v2/(Gamma__2*sigma__d), lambda__3 = Lambda__3*Gamma__1*sigma__v2/(Gamma__2*sigma__d3), sigma__v1 = Gamma__1*sigma__v2/Gamma__2}

Example - d(lambda__2)/d(sigma__d)

lambda2:=eval(eval(lambda__2,dims),Lambda__2=Lambda__2(Gamma__1,Gamma__2));

 

Lambda__2(Gamma__1, Gamma__2)*Gamma__1*sigma__v2/(Gamma__2*sigma__d)

Diff needs to know the functional dependencies of Gamma__1 and Gamma__2, but we withhold the actual function.

G1G2:={Gamma__1=Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2=Gamma__2(sigma__d, sigma__v2, gamma)};

{Gamma__1 = Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2 = Gamma__2(sigma__d, sigma__v2, gamma)}

Do the derivative.

deriv:=diff(eval(lambda2, G1G2),sigma__d);

((D[1](Lambda__2))(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*(diff(Gamma__1(sigma__d, sigma__v1, gamma), sigma__d))+(D[2](Lambda__2))(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*(diff(Gamma__2(sigma__d, sigma__v2, gamma), sigma__d)))*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d)+Lambda__2(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*(diff(Gamma__1(sigma__d, sigma__v1, gamma), sigma__d))*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d)-Lambda__2(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2*(diff(Gamma__2(sigma__d, sigma__v2, gamma), sigma__d))/(Gamma__2(sigma__d, sigma__v2, gamma)^2*sigma__d)-Lambda__2(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d^2)

Substitute the pieces we will do numerically with functions. Start with the D[1](Lambda__1) etc

DiGjsubs:={seq(seq(D[i](cat(Lambda__,j))(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))=fdiff(cat(Lnum,j),[i],[Gamma__1,Gamma__2]),j=1..2),i=1..2)}:
deriv:=subs(DiGjsubs,deriv);

(fdiff(Lnum2, [1], [Gamma__1, Gamma__2])*(diff(Gamma__1(sigma__d, sigma__v1, gamma), sigma__d))+fdiff(Lnum2, [2], [Gamma__1, Gamma__2])*(diff(Gamma__2(sigma__d, sigma__v2, gamma), sigma__d)))*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d)+Lambda__2(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*(diff(Gamma__1(sigma__d, sigma__v1, gamma), sigma__d))*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d)-Lambda__2(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2*(diff(Gamma__2(sigma__d, sigma__v2, gamma), sigma__d))/(Gamma__2(sigma__d, sigma__v2, gamma)^2*sigma__d)-Lambda__2(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d^2)

Next the functions Lambda__1 and Lambda__2 are replced by the numerical routines (in startup code).

deriv:=subs(Lambda__1(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))=Lnum1(Gamma__1,Gamma__2),
     Lambda__2(Gamma__1(sigma__d, sigma__v1, gamma), Gamma__2(sigma__d, sigma__v2, gamma))=Lnum2(Gamma__1,Gamma__2),deriv);

(fdiff(Lnum2, [1], [Gamma__1, Gamma__2])*(diff(Gamma__1(sigma__d, sigma__v1, gamma), sigma__d))+fdiff(Lnum2, [2], [Gamma__1, Gamma__2])*(diff(Gamma__2(sigma__d, sigma__v2, gamma), sigma__d)))*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d)+Lnum2(Gamma__1, Gamma__2)*(diff(Gamma__1(sigma__d, sigma__v1, gamma), sigma__d))*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d)-Lnum2(Gamma__1, Gamma__2)*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2*(diff(Gamma__2(sigma__d, sigma__v2, gamma), sigma__d))/(Gamma__2(sigma__d, sigma__v2, gamma)^2*sigma__d)-Lnum2(Gamma__1, Gamma__2)*Gamma__1(sigma__d, sigma__v1, gamma)*sigma__v2/(Gamma__2(sigma__d, sigma__v2, gamma)*sigma__d^2)

Now evaluate using the functional forms of Gamma__1 and Gamma__2

deriv:=eval(deriv,{Gamma__1(sigma__d, sigma__v1, gamma) = sigma__d*sigma__v1*gamma,
            Gamma__2(sigma__d, sigma__v2, gamma) = sigma__d*sigma__v2*gamma});

(fdiff(Lnum2, [1], [Gamma__1, Gamma__2])*sigma__v1*gamma+fdiff(Lnum2, [2], [Gamma__1, Gamma__2])*sigma__v2*gamma)*sigma__v1/sigma__d-Lnum2(Gamma__1, Gamma__2)*sigma__v1/sigma__d^2

Final scaling and removal of the gamma

scaledderiv:=eval(expand(sigma__d^2/sigma__v1*deriv), {sigma__v1=Gamma__1/sigma__d/gamma,sigma__v2=Gamma__2/sigma__d/gamma});

Gamma__1*fdiff(Lnum2, [1], [Gamma__1, Gamma__2])+fdiff(Lnum2, [2], [Gamma__1, Gamma__2])*Gamma__2-Lnum2(Gamma__1, Gamma__2)

Make it a function and check it out

scaled_dlambda_2dsigma_d:=unapply(scaledderiv,Gamma__1,Gamma__2);

proc (Gamma__1, Gamma__2) options operator, arrow; Gamma__1*fdiff(Lnum2, [1], [Gamma__1, Gamma__2])+fdiff(Lnum2, [2], [Gamma__1, Gamma__2])*Gamma__2-Lnum2(Gamma__1, Gamma__2) end proc

plot3d(scaled_dlambda_2dsigma_d(Gamma__1,Gamma__2),Gamma__1=0.1..5, Gamma__2=0.1..5,grid=[25,25]);

NULL

Download derivs2.mw

@mmcdara That's an interesting question. Since there is never a worksheet, I never respond (unless to comment after someone else has generated a worksheet.)

@acer Yes, this is the same idea. I originally was doing fsolve over the range, _Z = 0..infinity, which returns multiple real roots (but might miss some?), but I changed to fsolve( ,complex) to reliably return all roots and filtered out the complex and negative ones, and added fnormal and simplify(..,zero) as is my usual practice for complex roots.
 

Empirically Digits = 20 was an improvement, but I didn't try beyond that (the OP can try that). At least in the first case that still shows artifacts, the expression to be evaluated is much more complicated, so I'm guessing this is just propagation of errors. At Gamma__1=0, the expression is 0/0 so it is perhaps not surprising that near zero the numerical issues are more severe. I tried normal(,,,expanded) but it didn't make much difference. Probably some manipulation to a nicer form, e.g. continued fractions might work, but just brute force at higher Digits is less work.

@MaPal93 I wrote a preprocessing proc that automatically chooses the greatest RootOf - this has worked well for the first series of plots. For the second set there are still some artefacts (?) for small Gamma__1, which were reduced by going to Digits:=20; in preprocess. You could try higher Digits, but there might be some sort of singular behavior there.

derivatives_jumps_2.mw

@MaPal93 I just took the results from nondimoverall.mw above, which had T = sigma__v2/sigma__v1 and Gamma__1 = sigma__d*sigma__v1*gamma, and noticed that T*Gamma__1 could be called Gamma__2, replacing T and giving  Gamma__2 = sigma__d*sigma__v2*gamma and Gamma__1 = sigma__d*sigma__v1*gamma.

@MaPal93 I should have used the word greatest (most positive) rather than largest. Most positive real root meant the greatest positive one (we know there is at least one positive one). To prove uniqueness in a rigorous way is difficult, but just testing numerically over a grid shows that the solution exists everywhere and is unique.

nondimGamma1Gamma2Unique.mw

@MaPal93 

You tried two different non-dimensionalizations and concluded that in these two cases analytical/symbolic-form analysis leads to a Lambda_1 with some weird cliffs and a Lambda_2 that is NOT strictly positive, i.e., the surface has a boundary. In other words, you couldn't find a symbolic form for Lambda_1 and Lambda_2 that is consistent with the numerical solutions (at least for the first root of the degree-10 polynomial and for T and Sigma both not too close to 0).

These are not features of the symbolic solution, but of the numerical routines used by plot3d. If the RootOf for Lambda__1 is consistently interpreted as the most positive real root, and this value is used to calculate Lambda__2 then both Lambda__1 and Lambda__2 are consistently positive. The problem may be close roots not being properly resolved by the default numerical solver. Also, the order of the RootOfs if there is more than one positive root is to have the smallest one first.
The solution is to be careful about always choosing the most positive root for Lambda__1.

Questions :

1. (Assuming that it's tough to guess a non-dimensionalization Lambda_1 and Lambda_2 leading to two symbolic-forms whose plots are consistent with the well-behaved numerical solutions) is it possible to reverse-engineer a (even approximative) symbolic form from the numerical solutions?

See above. The existing symbolic solution is OK.

2. Even if we managed to reverse-engineer such symbolic form, would I encounter issues if I wanted to take partial derivatives of Lambda_1 and Lambda_2 wrt to sigma_v1 and sigma_v2 (since T and Sigma "entangles" both sigma_v1 and sigma_v2)? Specifically, you transitioned from a Gamma_1=sigma_v1*sigma_d*gamma and Gamma_2=sigma_v2*sigma_d*gamma non-dimensionalization to a T=sigma_v2/sigma_v1 and Sigma=gamma*sigma_d*sqrt(sigma_v1^2+sigma_v2^2) non-dimensionalization just to get rid of the square roots, but wouldn't be more "intuitive" to work with Gamma_1 and Gamma_2? (it's easier to interpret a 3D plot with Gamma_1 and Gamma_2 as x- and y- axes)

I reworked the solution non-dimensionalized in terms of Gamma__1 and Gamma__2, adding the better numerical evaluation. Since the dim -> non-dim and its reverse are just changes of variables, entangling should not be an issue - PDEtools:-dchange will change variables for partial derivatives.

3. Perhaps a stupid question: worst case scenario, are there ways to calculate partial derivatives when no analytical expressions are available? 

Yes, that's what fdiff does.

nondimGamma1Gamma2.mw

@Carl Love Aargh - thanks - I should have used indets.

@Steven_Huang Each problem is different in that respect. But you could just try all combinations - here there are 4 different solutions.

[Edit - fixed according to Carl's note]

restart

NULLN := (-t^2-1)*u^2-2*(t^2+1)^2*u+2*tNULL

(-t^2-1)*u^2-2*(t^2+1)^2*u+2*t

A := u^2*a[2]+u*a[1]+a[0]; B := u^2*b[2]+u*b[1]+b[0]

u^2*a[2]+u*a[1]+a[0]

u^2*b[2]+u*b[1]+b[0]

Q := collect(B*(diff(A, u))-A*(diff(B, u)), u)

(-a[1]*b[2]+a[2]*b[1])*u^2+(-2*a[0]*b[2]+2*a[2]*b[0])*u-a[0]*b[1]+b[0]*a[1]

vars := `minus`(indets(Q), {u})

{a[0], a[1], a[2], b[0], b[1], b[2]}

varcombs := combinat:-choose(vars, 3)

{{a[0], a[1], a[2]}, {a[0], a[1], b[0]}, {a[0], a[1], b[1]}, {a[0], a[1], b[2]}, {a[0], a[2], b[0]}, {a[0], a[2], b[1]}, {a[0], a[2], b[2]}, {a[0], b[0], b[1]}, {a[0], b[0], b[2]}, {a[0], b[1], b[2]}, {a[1], a[2], b[0]}, {a[1], a[2], b[1]}, {a[1], a[2], b[2]}, {a[1], b[0], b[1]}, {a[1], b[0], b[2]}, {a[1], b[1], b[2]}, {a[2], b[0], b[1]}, {a[2], b[0], b[2]}, {a[2], b[1], b[2]}, {b[0], b[1], b[2]}}

map2(solve, identity(N = Q, u), varcombs); nops(%)

{{a[0] = (t^4*a[1]+2*t^2*a[1]+2*t*a[2]+a[1])/(t^2+1), b[0] = (t^4*a[1]*b[1]+2*t^2*a[1]*b[1]+2*t^3+2*t*a[2]*b[1]+a[1]*b[1]+2*t)/((t^2+1)*a[1]), b[2] = (t^2+a[2]*b[1]+1)/a[1]}, {a[0] = (t^4*a[1]+2*t^2*a[1]+2*t*a[2]+a[1])/(t^2+1), b[0] = -(t^6-t^4*a[1]*b[2]+3*t^4-2*t^2*a[1]*b[2]-2*t*a[2]*b[2]+3*t^2-a[1]*b[2]+1)/((t^2+1)*a[2]), b[1] = -(t^2-a[1]*b[2]+1)/a[2]}, {a[0] = (t^4*a[1]+2*t^2*a[1]+2*t*a[2]+a[1])/(t^2+1), b[1] = -(-t^2*a[1]*b[0]+2*t^3-a[1]*b[0]+2*t)/(t^4*a[1]+2*t^2*a[1]+2*t*a[2]+a[1]), b[2] = (t^4+2*t^2+a[2]*b[0]+1)*(t^2+1)/(t^4*a[1]+2*t^2*a[1]+2*t*a[2]+a[1])}, {a[0] = (t^4+2*t^2+a[2]*b[0]+1)/b[2], a[1] = (t^6+3*t^4+t^2*a[2]*b[0]-2*t*a[2]*b[2]+3*t^2+a[2]*b[0]+1)/(b[2]*(t^4+2*t^2+1)), b[1] = (t^2*b[0]-2*t*b[2]+b[0])/(t^4+2*t^2+1)}, {a[0] = (t^4*a[1]*b[1]+2*t^2*a[1]*b[1]-2*t^3+2*t*a[1]*b[2]+a[1]*b[1]-2*t)/(b[1]*(t^2+1)), a[2] = -(t^2-a[1]*b[2]+1)/b[1], b[0] = (t^4*b[1]+2*t^2*b[1]+2*t*b[2]+b[1])/(t^2+1)}, {a[0] = (t^6+t^4*a[2]*b[1]+3*t^4+2*t^2*a[2]*b[1]+2*t*a[2]*b[2]+3*t^2+a[2]*b[1]+1)/(b[2]*(t^2+1)), a[1] = (t^2+a[2]*b[1]+1)/b[2], b[0] = (t^4*b[1]+2*t^2*b[1]+2*t*b[2]+b[1])/(t^2+1)}, {a[0] = -(-a[1]*b[0]+2*t)/b[1], a[2] = -(1/2)*(t^4*a[1]*b[1]-t^2*a[1]*b[0]+2*t^2*a[1]*b[1]+2*t^3-a[1]*b[0]+a[1]*b[1]+2*t)/(t*b[1]), b[2] = -(1/2)*(t^4*b[1]-t^2*b[0]+2*t^2*b[1]-b[0]+b[1])/t}, {a[0] = -2*t*(t^4+2*t^2+a[2]*b[0]+1)/(t^4*b[1]-t^2*b[0]+2*t^2*b[1]-b[0]+b[1]), a[1] = -2*(t^2+a[2]*b[1]+1)*t/(t^4*b[1]-t^2*b[0]+2*t^2*b[1]-b[0]+b[1]), b[2] = -(1/2)*(t^4*b[1]-t^2*b[0]+2*t^2*b[1]-b[0]+b[1])/t}, {a[0] = -(-a[1]*b[0]+2*t)*(t^4+2*t^2+1)/(t^2*b[0]-2*t*b[2]+b[0]), a[2] = -(t^6-t^4*a[1]*b[2]+3*t^4-2*t^2*a[1]*b[2]+3*t^2-a[1]*b[2]+1)/(t^2*b[0]-2*t*b[2]+b[0]), b[1] = (t^2*b[0]-2*t*b[2]+b[0])/(t^4+2*t^2+1)}, {a[1] = (a[0]*b[1]+2*t)/b[0], a[2] = -(1/2)*(t^4*a[0]*b[1]+2*t^5-t^2*a[0]*b[0]+2*t^2*a[0]*b[1]+4*t^3-a[0]*b[0]+a[0]*b[1]+2*t)/(t*b[0]), b[2] = -(1/2)*(t^4*b[1]-t^2*b[0]+2*t^2*b[1]-b[0]+b[1])/t}, {a[1] = (t^2*a[0]-2*t*a[2]+a[0])/(t^4+2*t^2+1), b[0] = (t^4*a[0]*b[1]+2*t^5+2*t^2*a[0]*b[1]+4*t^3+a[0]*b[1]+2*t)/(t^2*a[0]-2*t*a[2]+a[0]), b[2] = (t^2+a[2]*b[1]+1)*(t^4+2*t^2+1)/(t^2*a[0]-2*t*a[2]+a[0])}, {a[1] = (t^2*a[0]-2*t*a[2]+a[0])/(t^4+2*t^2+1), b[0] = -(t^4+2*t^2-a[0]*b[2]+1)/a[2], b[1] = -(t^6+3*t^4-t^2*a[0]*b[2]+2*t*a[2]*b[2]+3*t^2-a[0]*b[2]+1)/((t^4+2*t^2+1)*a[2])}, {a[1] = (t^2*a[0]-2*t*a[2]+a[0])/(t^4+2*t^2+1), b[1] = -(2*t^5-t^2*a[0]*b[0]+4*t^3+2*t*a[2]*b[0]-a[0]*b[0]+2*t)/(a[0]*(t^4+2*t^2+1)), b[2] = (t^4+2*t^2+a[2]*b[0]+1)/a[0]}, {a[1] = (a[0]*b[1]+2*t)*(t^2+1)/(t^4*b[1]+2*t^2*b[1]+2*t*b[2]+b[1]), a[2] = -(t^6+3*t^4-t^2*a[0]*b[2]+3*t^2-a[0]*b[2]+1)/(t^4*b[1]+2*t^2*b[1]+2*t*b[2]+b[1]), b[0] = (t^4*b[1]+2*t^2*b[1]+2*t*b[2]+b[1])/(t^2+1)}, {a[1] = (2*t^5+t^2*a[0]*b[0]+4*t^3-2*t*a[0]*b[2]+a[0]*b[0]+2*t)/((t^4+2*t^2+1)*b[0]), a[2] = -(t^4+2*t^2-a[0]*b[2]+1)/b[0], b[1] = (t^2*b[0]-2*t*b[2]+b[0])/(t^4+2*t^2+1)}, {a[2] = -(1/2)*(t^4*a[1]-t^2*a[0]+2*t^2*a[1]-a[0]+a[1])/t, b[0] = (a[0]*b[1]+2*t)/a[1], b[2] = -(1/2)*(t^4*a[1]*b[1]-t^2*a[0]*b[1]+2*t^2*a[1]*b[1]-2*t^3-a[0]*b[1]+a[1]*b[1]-2*t)/(t*a[1])}, {a[2] = -(1/2)*(t^4*a[1]-t^2*a[0]+2*t^2*a[1]-a[0]+a[1])/t, b[0] = 2*(t^4+2*t^2-a[0]*b[2]+1)*t/(t^4*a[1]-t^2*a[0]+2*t^2*a[1]-a[0]+a[1]), b[1] = 2*t*(t^2-a[1]*b[2]+1)/(t^4*a[1]-t^2*a[0]+2*t^2*a[1]-a[0]+a[1])}, {a[2] = -(1/2)*(t^4*a[1]-t^2*a[0]+2*t^2*a[1]-a[0]+a[1])/t, b[1] = -(-a[1]*b[0]+2*t)/a[0], b[2] = (1/2)*(-t^4*a[1]*b[0]+2*t^5+t^2*a[0]*b[0]-2*t^2*a[1]*b[0]+4*t^3+a[0]*b[0]-a[1]*b[0]+2*t)/(t*a[0])}}

18

 

Download Huang3.mw

@sursumCorda SCR submitted

Agree, something clearly very wrong here.

add(AllPairsDistance(G))/2

gives the correct result. (Note you need PathGraph(21) to get 20 edges.)

 

@sursumCorda @vv Thanks for the notes. I typically have done these directly with Minpoly and avoided MinimalPolynomial based on past experience though I can't remember why. I did the calculations here without MinimalPolynomial, and then added MinimalPolynomial as an afterthought. Looks like its result depends on the setting of Digits, which should not be the case IMO for what is supposed to be an exact calculation, which I guess is @sursumCorda's point.

@Joe Riel My mistake; I intended them to be the same size. (Usually I test and upload, but since Matrix output from Maple 2024 is not correctly rendered on Mapleprimes, I got lazy)

@goebeld The extend command was for older style matrices in the linalg package. These are no longer used (deprecated) as the help page tells you. Older matrices used "matrix", current ones use "Matrix" and commands in the LinearAlgebra package;

The brackets are entered with the < "less than" and > "greater than" keys, even though they display slightly differently. The objects produced are exactly the same as those from Matrix

The best way to extend a Matrix with zeros is probably just to embed it in a larger Matrix of zeroes.

XX := Matrix(2, 4, [[1, 1, 1, 1], [2, 2, 2, 2]]);
X2 := Matrix(5, 6); # Matrix of zeroes
X2[1..2,1..4] := XX:
X2;

 

First 24 25 26 27 28 29 30 Last Page 26 of 85