dharr

Dr. David Harrington

5550 Reputation

21 Badges

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

MaplePrimes Activity


These are replies submitted by dharr

@Nicole Sharp To apply fnormal and simplify/zero to the arguments of Quantity, you need a bit of advanced Maple. If q is the expression containing Quantity(..., ...),, then use

subsindets(q, specfunc(Quantity), x -> map(y->simplify(fnormal(y,23),zero), x));

This means every time we see the specific function Quantity, we apply the procedure given, which receives Quantity(a,b). The map over x does it for each of the ops of Quantity, which are the arguments a and b. For each of these, we apply the procedure y->simplify(fnormal(y,23),zero). For your case this yields

subsindets.mw

@C_R The numeric option for int, int(..., numeric), works identically to the eval/Int. (I think it was introduced to be a bit more intuitive than evalf/int and a similar syntax was already in place for dsolve.)  It happens that for the _Dexp method and these examples the results are all real.

Int_vs_numeric.mw

In one of Nicole's examples I tried all the high-digits methods and always got a complex result. This was for one with an infinite limit. The help page suggests there can be symbolic preprocessing for this case. There is also interaction with Units and Quantites-with-errors, which I think are more fickle (buggy?) than in their absence.

@Nicole Sharp I changed _s to Unit(s). Then numeric integration to infinity gives a complex value that fnormal can handle.

restart;

with(Units) :
Digits:=32;

Automatically loading the Units[Simple] subpackage
 

32

c := 299792458*Unit(m)/Unit(s) :
h := 662607015*10^(-8)*10^(-34)*Unit(J)/Unit(Hz) :
k := 1380649*10^(-6)*10^(-23)*Unit(J)/Unit(K) :
lambda_V_max := 750*Unit(nm) :
lambda_V_min := 380*Unit(nm) :
T_Sol := 5772.0*Unit(K) :
T0_Sol := 5772*Unit(K) :

evalf(Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T_Sol)) - 1)*lambda^5), lambda = 0*Unit(m) .. infinity*Unit(m),numeric = true));
fnormal(%);
simplify(%,zero);

(62938592.470335950467548474587301-0.31415926535897932384626433832795e-44*I)*Units:-Unit(kg/s^3)

(62938592.470335950467548474587301-0.*I)*Units:-Unit(kg/s^3)

62938592.470335950467548474587301*Units:-Unit(kg/s^3)

 

Download numeric_integration.mw

shoot expects the shooting parameters to be given in a list before output=listprocedure, something like [alpha = 0., beta =0.]. Actually I'm not sure if you can only have one parameter. I still didn't get it to work, and would need to understand the code better to figure that out. Perhaps author Doug Meade can help - he is on Mapleprimes

@Nicole Sharp I would expect, since the integrand doesn't involve any complicated functions, that the numeric integration always provides a real result. For me (version 2023.2.1) this is true even for the infinite limit: 

Digits:=32;
Pi*int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772)) - 1)*lambda^5), lambda = 0 .. infinity,numeric);

which gives

@dharr In this case, doing the integral numerically leads to a real answer differing in the last three places from evalf of the exact result at 32 digits

Digits:=32;
Mnum := (T, lambda_min, lambda_max) -> Pi*int(2*c^2*h/((exp(c*h/(k*lambda*T)) - 1)*lambda^5),
  lambda = lambda_min .. lambda_max, numeric) :
Mnum(T0_Sol, lambda_V_min, lambda_V_max);

@Nicole Sharp I agree; there certainly seems no reason that derive couldn't be any expression. And I also agree about updating the fundamental constants. Especially since many now have exact defined values that won't change in the future.

@Pemudahijrah01 ... and put a colon ":" at the end of the last line so you don't get your large output problem again. And use "local gamma;" at the beginning so gamma does not have a spcial meaning. (Like linalg, these were feedbacks you got on your earlier post.)

It seems that the form of the derive equation is very restrictive. From the help page, "In the derive_eqn equation of the form 'derive'=derive_obj, the derive_obj expression is typically a product of rational powers of numerics, Maple constants, and physical constant identifiers (for example, symbols). Exceptions are the abs function, and a sum with dimensionally consistent summands."
So perhaps you can only have one or the other and not a mix of these two types. If that is the reason, then there should at least be a warning or error message.

Here is a workaround:

restart;

with(Units[Standard]):with(ScientificConstants):

AddConstant(Solar_equatorial_radius, symbol = r[e,Sol], value = 696342., uncertainty = 65., units = km) ;

AddConstant(Solar_flattening, symbol = f[Sol], value = 0.000009) ;

AddConstant(Solar_eccentricity, symbol = e[Sol], derive = 1 - f[Sol]) ;

AddConstant(Solar_polar_radius, symbol = r[p,Sol], derive = r[e,Sol]*e[Sol]) ;

AddConstant(Solar_nonradius, symbol = x[Sol], derive = f[Sol]*r[e,Sol]) ;

evalf(Constant(r[e,Sol],units));

696342000.*Units:-Unit(m)

evalf(Constant(f[Sol],units));

0.9e-5

evalf(Constant(r[p,Sol],units));

696335732.9*Units:-Unit(m)

-696342e3*(0.000009-1)

696335732.9

evalf(Constant(x[Sol],units));

6267.078000*Units:-Unit(m)

NULL

Download constants.mw

@acer My understanding is that the minimal polynomial is the lowest degree polynomial (with rational coefficients since you didn't specify an extension field) that u1 with index=1 for both RootOfs is a root of. Since that u1 value is not the one corresponding to the desired trig expression, is it only coincidence that in this case the desired root is another root of the minimal polynomial? There is something special in this case that of the 6 roots only three are distinct, so does this work more generally? I suppose it can always be tried and checked...

@dharr I took a closer look at why @acer's solution doesn't lead to the solutions in the form with the square root, which Maple cannot simplify further.

restart;

expected:=(1+cos(arctan(3/4)/3))/3;evalf(%);

1/3+(1/3)*cos((1/3)*arctan(3/4))

.6590276223

alpha[1]:=evalc(convert(RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1,index=1),radical));evalf(%);
alpha[2]:=evalc(convert(RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1,index=2),radical));evalf(%);
alpha[3]:=evalc(convert(RootOf(60*_Z^3 - 60*_Z^2 + 15*_Z - 1,index=3),radical));evalf(%);

-(1/6)*cos((1/3)*arctan(3/4))+1/3-(1/6)*3^(1/2)*sin((1/3)*arctan(3/4))

.1090390090

-(1/6)*cos((1/3)*arctan(3/4))+1/3+(1/6)*3^(1/2)*sin((1/3)*arctan(3/4))

.2319333686

1/3+(1/3)*cos((1/3)*arctan(3/4))

.6590276223

The simplify here is the key to finding the simpler (factored) form of the quadratic). (a here stands in for the RootOf above).

quad[1]:=RootOf(simplify(eval(4*_Z^2 + (4*a - 4)*_Z + 4*a^2 - 4*a + 1,a=alpha[1])));
quad[2]:=RootOf(simplify(eval(4*_Z^2 + (4*a - 4)*_Z + 4*a^2 - 4*a + 1,a=alpha[2])));
quad[3]:=RootOf(simplify(eval(4*_Z^2 + (4*a - 4)*_Z + 4*a^2 - 4*a + 1,a=alpha[3])));

RootOf((cos((1/3)*arctan(3/4))-3*_Z+1)*(3^(1/2)*sin((1/3)*arctan(3/4))-cos((1/3)*arctan(3/4))-6*_Z+2))

RootOf(-(cos((1/3)*arctan(3/4))-3*_Z+1)*(3^(1/2)*sin((1/3)*arctan(3/4))+cos((1/3)*arctan(3/4))+6*_Z-2))

RootOf(36*_Z^2+(12*cos((1/3)*arctan(3/4))-24)*_Z+4*cos((1/3)*arctan(3/4))^2-4*cos((1/3)*arctan(3/4))+1)

The quadratic doesn't obviously factor in the "native" form, so using allvalues before substituting in the alpha's leads to pesky square roots.

factor(4*_Z^2 + (4*a - 4)*_Z + 4*a^2 - 4*a + 1);

4*_Z^2+4*_Z*a+4*a^2-4*_Z-4*a+1

Find all 6 roots

rts[1]:=[allvalues(quad[1])];evalf(%);
rts[2]:=[allvalues(quad[2])];evalf(%);
rts[3]:=[allvalues(quad[3])];evalf(%);

[-(1/6)*cos((1/3)*arctan(3/4))+1/3+(1/6)*3^(1/2)*sin((1/3)*arctan(3/4)), 1/3+(1/3)*cos((1/3)*arctan(3/4))]

[.2319333686, .6590276223]

[-(1/6)*cos((1/3)*arctan(3/4))+1/3-(1/6)*3^(1/2)*sin((1/3)*arctan(3/4)), 1/3+(1/3)*cos((1/3)*arctan(3/4))]

[.1090390090, .6590276223]

[-(1/6)*cos((1/3)*arctan(3/4))+1/3-(1/6)*(-3*cos((1/3)*arctan(3/4))^2+3)^(1/2), -(1/6)*cos((1/3)*arctan(3/4))+1/3+(1/6)*(-3*cos((1/3)*arctan(3/4))^2+3)^(1/2)]

[.1090390091, .2319333685]

rts[1][2];
rts[2][2];

1/3+(1/3)*cos((1/3)*arctan(3/4))

1/3+(1/3)*cos((1/3)*arctan(3/4))

 

Download RootOfX.mw

@nm Nice use of remove_RootOf; I didn't know about that. When I evalf([%]) your result I get 

[[0.6114019859 = 0., 0.2874388753 = 0., 0.1011591386 = 0.]],

but none of these correspond to the expected result, so I'm confused about that.

@fnavarro If I just evaluate the expression at Digits=40 I don't see a problem, there are no strange values in pts:

restart;
Digits:=40;
pts:=[seq([x,x*exp(-x)],x=0..20,0.2)]:
plots:-pointplot(pts,connect);

I agree if you look at the points in the plot structure for the regular plot without adaptive = true, they look strange, but as  @Preben Alsholm points out this is a function of the adaptive routine. Perhaps you have an specific example not involving plot.

@Art Kalb There are lots of non-commutative operators - "." and the ones starting with "&". I was thinking that since "." knew about inverses, b.b^(-1) =1, that it might be useful, but expand doesn't work (&* is special in this respect, since expand has some special code for &*). In the end I think you have to specify a group to make progress. If you are OK with specifying group elements as permutations, then that is the simplest way, since you can use most groups in the GroupTheory package, as in the following. (It would be nice if Elements(G) always had the identity first so I fiddled with it to force it.) If you want nicer symbols for the group elements you have to work a bit harder, as DrawCayleyTable does.

group-ring_2.mw

@Nicole Sharp If the Maple code is going to run on your local machine, then it needs to be downloaded, regardless of which file format is used. Import imports the whole contents of the file, and MapleCloud also uploads and downloads files. (Your Mapleprimes account should work on MapleCloud,  easiest to sign on from within Maple.) The files could be single functions, if you want to dowload just one of many functions. But I don't see how you can avoid re-downloading after an update.

The alternative would be to run the Maple code on a server, say one you have an account on. Then you are accessing always the current version of that file. Maplesoft also has a separate product Maplenet that allows users to run Maple code on a server through a web browser.

3 4 5 6 7 8 9 Last Page 5 of 56