dharr

Dr. David Harrington

8335 Reputation

22 Badges

21 years, 5 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 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;

 

@Scot Gould extend is in the old linalg package, so extend(XX,1,0,0) added 1 row and 0 cols of values 0, .i.e, a row of zeros.

@akbarhp As I said, the a(n,j):=coeff(eq[n],A[j]): line gives an error for n=6 and j=9 because coeff didn't work on eq[6], which contains an expression containing the word constcoeffsols, suggesting that constcoeffsols earlier in the code did not work. This is the same error message in the version you posted above. So you need to take a look at the printed out eq[6] and see if there is something being passed to constcoeffsols that doesn't look right. Or since the other ones worked, mabe this one is just too hard for constcoeffsols.

- if so you could just run the code up to n=5,j=5 and see if the code give sensible output for that smaller problem.

@akbarhp

Using more Digits than you need is poor practice, since it slows everything down, makes output hard to read, and you almost never need that accuracy. If you really, really, really need that accuracy, use it only after your code runs and you find you have accuracy issues.

Is n supposed to be an integer? Then use n:=2 and not n:=2.0. Note that you set X[n]: =  1.0, but then you use X[2], which is not 1.0. X[2] and X[2.0] are not the same. In general, avoid decimal points if you want an exact value. 

To create entries in a Matrix a, use a:=Matrix(16,16); before your loop, then assign the coefficients in the loop by

a[n,j]:=coeff(eq[n],A[j]):

then you just need  Alanda := simplify(a);

(I left this.)

If I add print(n,j) before the line 

a(n,j):=coeff(eq[n],A[j]):

I find if fails after n=6 and j=9.

So I now do

if n=6 and j=9 then print(n,j,eq[n],A[j]) end if;

then the long expression for eq[6] contains "constcoeffsols", meaning that that command did not work.

I suppressed the long output; commented out Digits:=100; and closed off the opening part of a procedure that wasn't completed in one execution group. At least in Maple 2024, it runs without any error messages. Perhaps you could provide a similarly shortened (without many pages of output) worksheet that shows the actual error message you describe.

Vib-code.mw

@acer I'm not sure I understand your point. The OP knows the matrix MM, having set it above in the code but didn't name it. So I'm assuming the definition of MM would be made before its first use. Presumably the OP knows eq13 is of the form v = MM %. v1 - MM %. v2, and wants it in the form v = MM %. simplify(v1-v2). 

@Christian Wolinski convert(..,Image) was introduced in Maple 2021.

First 26 27 28 29 30 31 32 Last Page 28 of 87