## dharr

Dr. David Harrington

## 6230 Reputation

19 years, 342 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

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

• The only real and positive solution is for lambda_2 = lambda_3, as shown explicitly in Nondimshort.mw. Is this solution the same as ans23[2] and ans32[2] in NondimshortAllEqns.mw (i.e., Lambda_2 = Lambda_3 = beta)?
Yes, if you look at the definition of beta in the alias, it is the same quartic as before.
• In Nondimshort.mw you type "Seems to always have one positive root - the curve Jumps at three places (these places can be found with solve(...,parametric).)". How to do so?
solve(Eqq111, Lambda__2, parametric);
Moreover, this curve is Lambda_2 as a function of positive Gamma...how to obtain the original lambda_2?
use the relationship lambda__2 = Lambda__2*sigma__v/sigma__d.
If I type solve(Eq2redox, lambda__2); map(simplify, [allvalues(%)]); I obtain 4 roots...which one is the positive one?
You cannot tell from looking at the solution, since it will depend on the values of the parameters. You have a better chance with solve(Eqq111, Lambda__2); map(simplify, [allvalues(%)]); because there is only one parameter Gamma, but even there it is very complicated to analyse. In principle, solve(Eqq111, Lambda__2, parametric,explicit); should give you the different solutions for different ranges of Gamma, but it only gives values at special values of Gamma, presumably because this case is too hard to solve. You can probably plot different solutions for the different ranges of Gamma between the critical values to figure this out.
Is it unique?
Descarte's rule of signs says there is one positive root, but which one it is will depend on the value of Gamma.

## remove units...

@Nicole Sharp Here is a worksheet with both your examples in which eval(q, Unit =1) works. Since you didn't provide the actual worksheet, it's hard to further diagnose. Or perhaps I have misunderstood your question?

RemoveUnits.mw

@dharr Here's the full case that I think better answers your question about the apparent asymmetry that arises from the way the equations are solved, but really the answers are not that different.

 >
 > local gamma;

# Consider a system of three equations:

 > eq1 := sigma__v^2/(lambda__1*(sigma__v^2/(2*lambda__1^2) + (5*sigma__d^2)/4)): eq2 := (gamma*sigma__v^2 + 2*lambda__3)*sigma__v^2/(((2*gamma*sigma__v^2 + 8*lambda__2)*lambda__3 + 2*gamma*lambda__2*sigma__v^2)*(2*(gamma*sigma__v^2 + 2*lambda__3)^2*sigma__v^2/((2*gamma*sigma__v^2 + 8*lambda__2)*lambda__3 + 2*gamma*lambda__2*sigma__v^2)^2 + lambda__2^2*(gamma*sigma__v^2 + 4*lambda__3)^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + gamma^2*sigma__v^4*lambda__3^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__2)*lambda__3 + 2*gamma*lambda__2*sigma__v^2)^2 + sigma__d^2)): eq3 := (gamma*sigma__v^2 + 2*lambda__2)*sigma__v^2/(((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)*(2*(gamma*sigma__v^2 + 2*lambda__2)^2*sigma__v^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + lambda__3^2*(gamma*sigma__v^2 + 4*lambda__2)^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + gamma^2*sigma__v^4*lambda__2^2*sigma__d^2/((2*gamma*sigma__v^2 + 8*lambda__3)*lambda__2 + 2*gamma*lambda__3*sigma__v^2)^2 + sigma__d^2)):
 > Eq1, Eq2, Eq3 := eq1 - lambda__1, eq2 - lambda__2, eq3 - lambda__3:

"non-dimensionalize" Eq 1 by hand. New non-dim variable Lambda__1. End up with no parameters in the equation

 > Lambda__1 = lambda__1*sigma__d/sigma__v; isolate(%,lambda__1); Eqq1:=simplify(eval(Eq1,%))*sigma__d/sigma__v; solve(Eqq1,Lambda__1);

We can find the non-dimensionalization here by a systematic process, but let's guess Lambda__2 and Lambda__3 are  defined similarly to Lambda__1

 > {Lambda__2 = lambda__2*sigma__d/sigma__v,Lambda__3 = lambda__3*sigma__d/sigma__v}; ss:=solve(%,{lambda__2,lambda__3}); Eqq2:=numer(normal(eval(Eq2,ss),expanded)); Eqq3:=numer(normal(eval(Eq3,ss),expanded));

We can divide through by sigma__v, which is in every term, and we see that Gamma = gamma*sigma__v*sigma__d

 > Gamma = gamma*sigma__d*sigma__v; sg:=isolate(%,gamma); Eqqq2:=simplify(eval(Eqq2/sigma__v,sg)); Eqqq3:=simplify(eval(Eqq3/sigma__v,sg));

 > alias(alpha=RootOf((Gamma^2 + 4)*_Z^2 + 4*Gamma*_Z + Gamma^2), beta=RootOf(40*_Z^4 + 36*Gamma*_Z^3 + (9*Gamma^2 - 4)*_Z^2 - 4*Gamma*_Z - Gamma^2), delta=RootOf((200*Gamma^2 + 800)*_Z^6 + 560*Gamma*_Z^5 + (-24*Gamma^2 - 80)*_Z^4 + (-52*Gamma^3 - 64*Gamma)*_Z^3 + 16*Gamma^2*_Z^2 + 24*Gamma^3*_Z + 5*Gamma^4));

Solve in both orders - we see that we have the Lambda__2=Lambda__3 solution that we saw before. Lambda__2 and Lambda__3 are swapped, i.e. the apparent complexity depends on the solving order

Consider the third solution where one solution looks simpler than the other. But after taking allvalues  we see that finally the complexity is the same although it didn't look like it. For this third solution the Lambda's are different (complex conjugates) and complex

 > ans23[3]; map(simplify,[allvalues(%)]); ans32[3]; map(simplify,[allvalues(%)]);

 >

## dimensional analysis...

@dharr For the case where lambda__2 = lambda__3, there is always a positive root. The same analysis could also be done for the case where the two lambda's are not set equal.

Edit - actually Descartes rule of signs tells us there must be one positive root no matter what the sign of the term with coefficient 4-9*Gamma^2

Nondimshort.mw

## Coments...

2. In general I understand, but I thought my commands l1-remove(has, l1, I), l2-remove(has, l2, I), and l3-remove(has, l3, I) would return the complex solutions for the respective lambda, if any exist. Since they all return zero, I thought this meant all solutions were real. Am I wrong?

has(...,I) is not a good test. Firstly, an expression may have some I's that cancel if the expression is simplify using evalc. more importantly, everything depends on the parameters. Sor example sqrt(1-sigma__v) might be real or complex depending on the values of sigma__v even though there are no explicit I's.

3. I am not sure I got you. In your example the two specular equations generate specular elements for a solution. Don't look at your av variable but at your ans variable (i.e., implicit form) instead. In particular, look at your third one (with RootOf):

They are not exactly the same, but the RootOf term is the same (a quadratic). In short, the difference is just a sign and the -1. Instead, I get:

Check output 16 of the new script I already attached to see more clearly.

It is not the convolutedness that surprises me, but the fact that lambda_3 is a quartic in _Z, while lambda_2 is made of some weird combination of (actually the same) quartic and with some randomly looking leading coefficients. My question: is this "asymmetry" plausible from two perfectly symmetric equations? In particular, how can the significant difference in the size of the implicit forms for lambda_2 and lambda_3 be explained?

Respone: I'm not sure which file has the relevant eq 16 in; not the last one you posted?. I expect that if you have the quadratic formula of the  quartic formula, it might look different from the quartic formula on the quadratic formula, even if they are actually the same. (I would solve the quadratic first becuase is is easier to see what the right root might be.) To see if they are the same, choose some random values of the parameters and see if they actually evaluate to the same number.

## removing units...

@Nicole Sharp There are definitely problems with the Units package that Maple is working on. In your last two cases though, the Units package just sees a Maple function Quantity() that if specified might return something with units that weren't the same as those of its arguments.

If q is an expression with units and you want to remove them all, just use eval(q, Unit =1)

1. What am I doing wrong when solving the reduced form system?
Nothing, the expressions are just too complicated. You substitute the solution of a quartic into a complicated eqn and then that eq is very complicated.
2. Am I understanding correctly that, in contrast to the provided example in the answer, only real solutions exist here (even if lambda_2 \= lambda_3)?
3. lambda_2 (l2[1]) and lambda_3 (l3[1]) are pretty much nonspecular (i.e., very different in size and form), while in the example provided in the answer the solutions are specular/very similar to each other, am I wrong?
In my solution, the stepwise equation get more complicated (admittedly only slightly). I think this is general. Of course if you have the same values they potentially can be simplified to the same thing, but once there are radicals, then that is not easy to do. You also don't know that you have picked out the right roots; they might be different.

You are assuming the first value from allvalues is real, but this will depend on your parameters. Only for numeric coefficients is the first RootOf real (if there is a real root).

240124_specular_equations_nonspecular_solutions_dharr.mw

## reason...

Presumably, the reason that f(0.2342493224)  is not changed to f(0.2342)  is because to get 4 place accuracy in the value of the function might require more than 4 place accuracy in its argument.

## seq...

@acer Thanks. Obvious in highsight, I should have gone straight to seq.

I think that getting different results depending in units present or not does seem strange, as does different results for numbers being expressed in different forms (exact or float). These differences seem to be becuase there is different analysis going on behind the scenes. Those can certainly be considered areas for improvement, but I don't think something can be considered a bug unless it returns the wrong answer. A numerical calculation is expected to have a small error. Since Maple works in the complex domain a small imaginary part is no different from a small error in a real answer, which is why the fnormal and simplify/zero postprocessing seems to me to be the correct way to go.

A similar comment applies to "There should be an option in the Maple settings to instruct Maple that all integrations should only give real results and thus Maple can then automatically try alternative integration methods whenever a non-real result is detected." But what you really mean is "if the integrand is always real then the integral should return real", or otherwise you do not allow for integration involving complex quantities. But a small change to the arguments of some functions in the integrand can make them complex, and then you are back to the same problem. The RealDomain package tries to solve this for some functions (but not int), and does not seem to work as well as might be expected.

There is also the possibility of trying different methods to get a real result. But for a really difficult problem, the time taken for this might be quite long, so you would want this only as an option. Many Maple routines do allow you to specify a list of methods or exclude some methods or exclude symbolic preprocessing (which would perhaps help most here). But for integration, the symbolic part can be powerful, say in detecting singularities that might otherwise be missed and that might be degrading accuracy.

## fnormal inside Quantity...

@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

## Int vs int(..., numeric)...

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

## numeric integration...

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

 > 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);

 >

## syntax...

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

## numeric integration...

@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

 First 10 11 12 13 14 15 16 Last Page 12 of 64
﻿