dharr

Dr. David Harrington

8832 Reputation

22 Badges

21 years, 121 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

@Aung If you ask FunctionAdvisor this

FunctionAdvisor(special_values, hypergeom([1/2, -k/2], [1 - k/2], z))

and look through the possibilities, most require special values such as integer of half-integer for the parameters. The exception seems to be for JacobiP, and indeed

convert(y, JacobiP) assuming k::positive;

gives

but I'm guessing that is not what you want. I agree with @Preben Alsholm that there is no reason not to use the hypergeom function: it can be evaluated, plotted, differentiated etc just like any other function.

@MaPal93 Here's some progress, though Lambda__2 is so complicated it is very slow to find the boundary between positive and negative values, e.g, by implicitplot or even fsolve. I'll give some thought to a faster method. I didn't explore the other solutions to see if the first RootOf is the only one.

Overall.mw

Edit: In principle, substituting Lambda__2 =0 into the original equations and solving is simpler than working with the complicated solutions, but doesn't seem to agree - maybe it is implicitly using the wrong solutions. I also did the non-dim in one step.

Overall2.mw

Its possible there is some better non-dimensionalization that is simpler to work wth, but I guess somehow the degree 10 polynomial is essential.

Edit2: Here's a different non-dimensionalization but the conclusions are the same.

nondimoverall.mw

@mmcdara I'm not sure what was happening there; it runs as expected in 2024.

@Ronan Sorry, I've uploaded the correct file now.

Here's some variations on @Scot Gould's answer. In general, use sort to solve these session-to-session variations.

MaplePrimes_sort_eigenvectors_and_eigenvalues.mw

 

@Ronan I don't use initialization files so you'll have to look at the documentation for where to find them. But if the output comes in a fresh session after restart but not until infolevel[RationalTrigonometry]:=0; it suggests that the mention of RationalTrigonometry is somehow loading it from the .mla file. userinfo allows you to use any name, so try userinfo(2, mytest,...) and infolevel[mytest]:=3;

@Ronan Your output suggests the ModuleLoad occurs before the with() statement? Are you loading it in some initialization file? I copied your userinfo() to my module and it works fine, so no syntax errors.

Otherwise, I'm out sof suggestions, sorry.

@Ronan userinfo is an executable statement, so should be after the uses and global declarations of the ModuleLoad procedure. Can you confirm that ModuleLoad is actually running by putting a print statement in? 

@Carl Love But from the user's (OP's?) point of view, if they have row and column vectors already, and want to make the rank 1 matrix, then c.r is the natural way to do it, which is consistent with the mathematical use of the `.` operation for vectors and matrices. The fact that it might also be an outer product and Maple has programmed it that way should be hidden when using `.`. 

@dharr Since the ModuleLoad: in the message may not be what you want, you can rename it to something useful, e.g.,

pkg := module() export squareme; option package, load = `Default global settings`; local `Default global settings`;
  `Default global settings` := proc() userinfo(2,pkg, "All things are blue"); end proc;
  squareme:=proc(x) x^2 end proc;
end module; 

 

@dharr It turns out if you set the initial conditions to be arbitrary functions of beta__c that reduce to the classical values at beta__c = 0, then you can find what is probably the most general solution.

restart;

ode1 := diff(c(T),T)=-2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2);
ode2 := diff(p__c(T),T)=2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T);
sys := {ode1,ode2}:

diff(c(T), T) = -2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2)

diff(p__c(T), T) = 2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T)

We want a solution that tends to the following solution as beta__c tends to zero.

sys0:=eval(sys,beta__c=0);
ans0:=dsolve(sys0 union {p__c(0)=pc0,c(0)=c0});

{diff(c(T), T) = -2*c(T), diff(p__c(T), T) = 2*p__c(T)}

{c(T) = c0*exp(-2*T), p__c(T) = pc0*exp(2*T)}

Let's define f(beta__c) as any function that has f(0) =c0 and g(beta__c) as any function that has g(0) = pc0.

f(0):=c0; g(0):= pc0; # remember this

c0

pc0

 

Solve - there are many solutions depending on signs of various things. Probably could work through the correct case(s),

but here I just use the first solution and make assumptions that make it look nice.

ans1:=dsolve(sys union {p__c(0)=g(beta__c),c(0)=f(beta__c)},[p__c(T),c(T)]):
nops([ans1]);
ans11:=simplify(ans1[1]) assuming beta__c>0, f(beta__c)>0, g(beta__c)>0, T>0;
simplify(eval(normal(%),beta__c=0)) assuming f(beta__c)>0, g(beta__c)>0, T>0, pc0>0:
ans11beta0:=simplify(%) assuming pc0>0;

8

{c(T) = f(beta__c)*g(beta__c)^(1/2)/((beta__c*f(beta__c)^2+g(beta__c)^2)*exp(8*T)-beta__c*f(beta__c)^2)^(1/4), p__c(T) = ((beta__c*f(beta__c)^2+g(beta__c)^2)*exp(8*T)-beta__c*f(beta__c)^2)^(1/4)*g(beta__c)^(1/2)}

{c(T) = c0*exp(-2*T), p__c(T) = pc0*exp(2*T)}

odetest(ans11,sys);

{0}

Try the other solving order - solution 8 is the same

ans2:=dsolve(sys union {p__c(0)=g(beta__c),c(0)=f(beta__c)},[c(T),p__c(T)]):
nops([ans]);
ans28:=simplify(ans2[8]) assuming beta__c>0, f(beta__c)>0, g(beta__c)>0, T>0;
simplify(eval(normal(%),beta__c=0)) assuming beta__c>0, f(beta__c)>0, g(beta__c)>0, T>0:
simplify(%) assuming pc0>0

1

{c(T) = f(beta__c)*g(beta__c)^(1/2)/((beta__c*f(beta__c)^2+g(beta__c)^2)*exp(8*T)-beta__c*f(beta__c)^2)^(1/4), p__c(T) = ((beta__c*f(beta__c)^2+g(beta__c)^2)*exp(8*T)-beta__c*f(beta__c)^2)^(1/4)*g(beta__c)^(1/2)}

{c(T) = c0*exp(-2*T), p__c(T) = pc0*exp(2*T)}

odetest(ans28,sys);

{0}

NULL

Download dsolve_simple.mw

 

@srmaple I agree that it seems Maple should have found a solutiion that changed smoothly as beta__c tended to zero. On the other hand it seems that you were unlucky in some ways that might not happen with other examples. When Maple eliminated the first variable (second in order given) it gave an equation without beta_c for the second variable, and that happened for both solving orders. Then the solution it did find had a beta__c in the denominator which meant dealing with a 0/0 form.

@srmaple Actually, I was comparing with the beta__c=0 solution, not Mathematica's, and suggested c2 = c2*beta__c as a simple function that went to zero. I even suggested c2*beta__c^2 was OK. But perhaps this is a better justification?

restart;

with(DEtools):

ode1 := diff(c(T),T)=-2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2);
ode2 := diff(p__c(T),T)=2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T);
sys := [ode1,ode2]:

diff(c(T), T) = -2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2)

diff(p__c(T), T) = 2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T)

First solve the system with beta__c = 0.
We want the solution for non-zero beta__c to reduce to this in the limit of beta__c -> 0.

sys0:=eval(sys,beta__c=0);
dsolve(sys0);

[diff(c(T), T) = -2*c(T), diff(p__c(T), T) = 2*p__c(T)]

{c(T) = c__2*exp(-2*T), p__c(T) = c__1*exp(2*T)}

sol2:=simplify([eval(dsolve(sys,[c(T),p__c(T)],explicit),c__2=f(beta__c))]) assuming beta__c>0;

[{c(T) = -2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = -2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = -2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = -2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}, {c(T) = 2*2^(1/2)*(-f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*f(beta__c))^(1/4)}]

So p__c(T) already tends to the right solution for f(beta__c) -> 0.

What other condition on f(beta) is required for c(T)? Take the first solution for example

cT1:=eval(c(T),sol2[1]);

-2*2^(1/2)*(f(beta__c)/(2*c__1*exp(8*T)+16*f(beta__c))^(1/2))^(1/2)/beta__c^(1/2)

Take a Taylor series for f(beta__c) with f(0)=0

ff:=eval(series(f(beta__c),beta__c),f(0)=0);

 

series((D(f))(0)*beta__c+((1/2)*((D@@2)(f))(0))*beta__c^2+((1/6)*((D@@3)(f))(0))*beta__c^3+((1/24)*((D@@4)(f))(0))*beta__c^4+((1/120)*((D@@5)(f))(0))*beta__c^5+O(beta__c^6),beta__c,6)

series(eval(cT1,f(beta__c)=ff),beta__c,2) assuming beta__c::positive;
eval(convert(%,polynom),beta__c=0); #first term

series(-2*2^(1/4)*((D(f))(0)/(c__1*exp(8*T))^(1/2))^(1/2)+((1/2)*2^(1/4)*((D(f))(0)/(c__1*exp(8*T))^(1/2))^(1/2)*(-((D@@2)(f))(0)*c__1*exp(8*T)+8*(D(f))(0)^2)/(c__1*exp(8*T)*(D(f))(0)))*beta__c+O(beta__c^2),beta__c,2)

-2*2^(1/4)*((D(f))(0)/(c__1*exp(8*T))^(1/2))^(1/2)

Above agrees with the beta__c=0 solution for any function of beta__c that goes to zero and whose first derivative is non-zero at beta__c=0, e.g., const*beta__c

NULL

Download dsolve.mw

@Ronan A couple of comments

1. I wasn't using makehelp, but see here for how to do subfolders with that method

Creating Help - An Exploration Of Makehelp The browser entries have 3 items. But notice that he did not succeed in altering the order.

2. I was supposing that you might just replace the TOC with a new one, by making the appropriate Record structure. In my case I wanted my Orbitals package under Physics, so that the the top level Record is for Physics, and it has a child Record of Orbitals, which has child Records for the individual help pages. Here's what the final TOC looks like, from lprint(eval(TOC)):

Record('entry' = Physics,'topic' = (NULL),'priority' = (NULL),'language' = "en"
,'product' = "Maple",'category' = "Help Page",'children' = [Record('entry' = 
Orbitals,'topic' = (NULL),'priority' = (NULL),'language' = "en",'product' = 
"Maple",'category' = "Help Page",'children' = [Record('entry' = "Overview",'
topic' = "Orbitals,Overview",'priority' = 100,'language' = "en",'product' = 
"Maple",'category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(atom4e, atomJ, atomK) - atomic exchange, coulomb and four-electron integrals"
,'topic' = "Orbitals,atom4e",'priority' = 70,'language' = "en",'product' = 
"Maple",'category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(cartesian, fullcartesian) - convert to cartesian coordinates",'topic' = 
"Orbitals,cartesian",'priority' = 65,'language' = "en",'product' = "Maple",'
category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(HYD) - hydrogenic orbital",'topic' = "Orbitals,HYD",'priority' = 60,'language
' = "en",'product' = "Maple",'category' = "Help Page",'children' = (NULL)), 
Record('entry' = "(HYDr) - radial part of hydrogenic orbital",'topic' = 
"Orbitals,HYDr",'priority' = 55,'language' = "en",'product' = "Maple",'category
' = "Help Page",'children' = (NULL)), Record('entry' = 
"(overlap) - overlap integral",'topic' = "Orbitals,overlap",'priority' = 50,'
language' = "en",'product' = "Maple",'category' = "Help Page",'children' = (
NULL)), Record('entry' = "Plotting orbitals examples",'topic' = 
"Orbitals,plots",'priority' = 5,'language' = "en",'product' = "Maple",'category
' = "Help Page",'children' = (NULL)), Record('entry' = 
"(realY) - real combinations of spherical harmonics",'topic' = "Orbitals,realY"
,'priority' = 45,'language' = "en",'product' = "Maple",'category' = "Help Page"
,'children' = (NULL)), Record('entry' = "Orbital package references",'topic' =
"Orbitals,references",'priority' = 2,'language' = "en",'product' = "Maple",'
category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(STO) - Slater type orbital",'topic' = "Orbitals,STO",'priority' = 35,'
language' = "en",'product' = "Maple",'category' = "Help Page",'children' = (
NULL)), Record('entry' = "(STOr) - radial part of Slater type orbital",'topic'
= "Orbitals,STOr",'priority' = 30,'language' = "en",'product' = "Maple",'
category' = "Help Page",'children' = (NULL)), Record('entry' = 
"(Y) - spherical harmonic",'topic' = "Orbitals,Y",'priority' = 25,'language' =
"en",'product' = "Maple",'category' = "Help Page",'children' = (NULL)), Record(
'entry' = "Orbitals package examples",'topic' = "Orbitals,examples",'priority'
= 10,'language' = "en",'product' = "Maple",'category' = "Help Page",'children'
= (NULL))])])

If you think it's useful, I could send the help files themselves, so you could try converting them to the help database.

@MaPal93 Your treatment neglects the fact that the relationship between your new Sigma_v and sigma__v1 and sigma__v2 is lost if you only proceed with p3 and don't add the new auxilliary equation to the analysis. I fixed this and the scale factor up.

For two equations convert each to polynomials by removing the radicals. Then each polynomial is normalized so the first terms is 1 and then all the remaining terms for the two polynomials and the auxilliary equation(s)  are collected (even if there is a lot of redundancy, the matrix procedure will take care of it.)

Case_with_radicals.mw

First 35 36 37 38 39 40 41 Last Page 37 of 95