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

@Carl Love I misinterpreted that statement as meaning you could have two different canonical orders even if the graphs were isomorphic. But I see your interpretation is correct. But the description is also correct; op(4,G) are the same for those graphs; the same problem also exists with @lcz's original list. The error has to do with a misplaced op(1,..) and perhaps 'nolist' for entries (always a problem!). One solution is 

op~(1,{entries}(ListTools:-Classify(
    g-> [seq](op(4, GraphTheory:-CanonicalGraph(g))),graph_list),'nolist'));

 

@Carl Love Yes, I use it when I can find such a procedure. In this case, it needs some work; the example is from the CanonicalGraph help page.

Download classify.mw

@lcz To add to @Carl Love 's explanation, op(4,G) is the same information you would get from Neighbors(G), but with generic labelling, and an Array of sets rather than a list of lists. See ?graph,details.

@mmcdara Very nice. Vote up.
@KIRAN SAJJAN Note that the paper has -Gr/(R*Re)*H' and not +Gr/(R*Re)*H'

@mmcdara Thanks! I was still bothered by not taking advantage of duplicate rows, and have now uploaded  above a version that can manage the L=10 case in about 1 s on my machine.

Since there is more likely not a match and that can be detected at the first mismatch, you want to try to not to generate all the sequence unless you have to (as in the McCarthy rules for 'and'). I don't think the timing tests are that reliable here. I suspect there are greater efficiencies to be found in the larger algorithm, so this is just a reply, not an answer.

tests.mw

@janhardo The error message had "unexpected option: [-50 .. 50, -5 .. 5, 0 .. 5] = [-50 .. 50, -5 .. 5, 0 .. 5]", where the option should have been "view = [-50 .. 50, -5 .. 5, 0 .. 5]", which suggests that the lhs view had unexpectedly obtained a value. That in turn suggests that the view parameter for complex_surface should be named something else. Another fix would be to have the complexplot3d option as ':-view' = view but I think renaming to something else avoids confusion.

@janhardo I fixed some small things. The error was fixed by using a different name than view for the complex_surface parameter. You weren't returning the generated plot because of the following printf statement (added to debug?).

restart;

complex_surface := proc(f, z_range::range(complex), plotview::list) local plt;
    
    printf("Complex function f(z) defined.\n");
    printf("Range of z values: %a\n", z_range);
    
    # Plot de complexe functie met complexplot3d
    printf("Plotting complex function...\n");
    plt := plots:-complexplot3d(f, z_range, 'view' = plotview, 'grid' = [50, 50]);
    
    printf("Plot generated.\n");
    
    plt   # return result
end proc:

# Definieer de complexe functie f(z) = ln(z)
f := z -> ln(z);

# Definieer het bereik van z-waarden voor de plot
z_range := -50-5*I .. 50+5*I;

# Roep de procedure aan om de plot te genereren
complex_surface(f, z_range, [-50..50, -5..5, 0..5]);

proc (z) options operator, arrow; ln(z) end proc

-50-5*I .. 50+5*I

Complex function f(z) defined.
Range of z values: -50-5*I .. 50+5*I
Plotting complex function...
Plot generated.

NULL

NULL

NULL

Download complexplot.mw

You might consider complexplot3d, e.g., 

plots:-complexplot3d(Zeta(z), z=-50-5*I..50+5*I,view=[-50..50, -5..5, 0..5],
                    grid=[50,50]);

gives

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

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