Carl Love

Carl Love

28035 Reputation

25 Badges

12 years, 320 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Form-conversion commands often need to be composed. The evalc command will take most expressions and express them as the sum of their real and imaginary parts. That is, the expression is put into cartesian or rectangular form. When polar is passed a single expression in this form, it converts it to polar form. Then evalf converts to decimal form.

polar(evalc(z));

     polar(sqrt(2), Pi/4)

evalf[4](%);

     polar(1.414, 0.7854)

The command randpoly might be a place to start:

randpoly(x);
randpoly(x, degree= 2);

etc.

Use literally the word undefined. Maple knows what this means.

A simpler way to solve the problem is to pass solve an inequality that forces the solution to be positive. This avoids the need for select, is, or another command line.

f:= gamma*sqrt(4) * sqrt(omega^2*C2^2*R4^2 /
     (C2^4*R4^4*beta^2*gamma^2*omega^4+C2^2*(1+gamma^2*(beta+1)^2-2*gamma) *
      R4^2*omega^2+1)):

solve({diff(f,omega), omega > 0}, omega) assuming positive;

     {omega = 1/(C2*R4*sqrt(beta*gamma))}

I have a hunch that the above is more robust than the flaky and ancient is because it uses the more modern and more sophisticated package SolveTools:-Inequality. However, I also know that it's still easy to confuse solve with inequalities, especially when square roots are involved.

An equivalent command is

solve(diff(f,omega), omega, useassumptions) assuming positive;

It's important for you to understand that the roots that you originally computed with RootFinding:-Isolate (or, equivalently, fsolve) are accurate to 10 digits. The large errors are introduced by your attempt to compute the residuals. Kitonum's Answer may give you the false impression that you need to increase Digits to compute the roots, although that wasn't his intention. No, you only need to increase Digits for the computation of the roots if you also want to compute the residuals, which isn't a necessary part of the computation of the roots.

restart:

p:=
     -12116320194738194778134937600000000*t^26 +
     167589596741213731838990745600000000*t^24 +
     1058345691529498270472972795904000000*t^22 -
     4276605572538658673086219419648000000*t^20 -
     23240154739806540070988490473472000000*t^18 -
     5442849111209103187871341215744000000*t^16 +
     49009931453396028716875310432256000000*t^14 +
     74247033158233643322704589225984000000*t^12 -
     2762178990802317464801412907008000000*t^10 -
     25947900993773120244883450232832000000*t^8 -
     7468990043547273070742668836864000000*t^6 -
     567730116675454293925108383744000000*t^4 +
     3703566799705707258760396800000000*t^2 -
     4742330812072533924249600000000
:

Some exact computations can be used to simplify your polynomial. These aren't necessary.

p1:= p/icontent(p):   #Remove common integer factor.

p2:= algsubs(t^2= T, p1):   #The polynomial is even.

p3:= unapply(convert(p2, horner), T)@(t-> t^2):   #Reduce amount of computation for residuals.

The rest of the computations are floating point.

Digits:= 10:

display:= proc(e, digits::posint) print(evalf[digits](e)); e end proc:

Compute the positive roots.

R10:= display([fsolve(p2, T= 0..infinity)^~(1/2)], 10):

[0.4212510777e-1, 0.6596008501e-1, .8048430445, 1.300314688, 2.295186769, 4.162501845]

Compute the residuals.

display(p3~(R10), 2):

[0.2e8, 0.1e8, -0.67e13, -0.32e17, 0.14e22, 0.11e27]

Now we do the exact same computations at 50 digits.

Digits:= 50:

R50:= display([fsolve(p2, T= 0..infinity)^~(1/2)], 10):

[0.4212510777e-1, 0.6596008501e-1, .8048430445, 1.300314688, 2.295186769, 4.162501845]

Note that the above roots are exactly the same (to 10 digits) as the previously computed roots.

display(p3~(R50), 2):

[0., 0.3e-32, -0.12e-26, 0.32e-23, -0.92e-19, 0.42e-13]

...but the residuals computed at 50 digits are much smaller in magnitude. It's also interesting that as we increased Digits by 40, the magnitudes were reduced by a factor of about "10^(40)." That's typical, but not guaranteed.



Download accuracy_of_roots.mw

To say that a matrix is "augmented", in this context, means that the right-side vector b has been appended as the rightmost column of the matrix. To achieve this, you simply need to make the command

IterativeApproximate(<A|b>, ...);

So, the b is incorporated into the first argument. In this case, don't use b as the second argument.

I get this solution:

     Vector[column](4, [3.98912386882810, 2.99292198093478, 2.99317583909596, 2.99337057407024])

which has residual

LinearAlgebra:-Norm(A.% - b);

     0.103910640752635786e-2

which is very close to your specified tolerance of 1e-3.

The above is a response to your originally posted Question, using b = <1,1,1,1>.

Responding to your reformulated Question, using b = <1,0,0,0>: The solution that you got seems correct to me. The residual is 0.235228835640644007e-3, which is well less than your specified tolerance.

 

It's important for a beginner to note that print doesn't "give" anything that can be stored for later use. It's often difficult for a beginner to understand this. The print command only displays something on the screen; it's rarely used in formal Maple programming.

Here's how to create A in one line of code and then create in a second line of code:

Download Matrices.mw
A:= Matrix(

     [[0$7], [0$6], [a,-d,a,-d,a], [-a,d,d,a], -[d$3], [d,d]],
     scan= band[0,5],
     shape= antisymmetric
);

C:= [seq(seq(eval(A, {a= i, d= j}), i= [-1,1]), j= [-1,1])];

Now C is a list of four matrices. The matrices can be accessed as C[1]C[2], etc. Note that for loops can often be replaced by seq expressions, and it's usually beneficial to do so.

Example:

A:= Array([4,3,7]);

If your objective function is a procedure, then the constraints must also be procedures. Your current constraints are expressions. So your command should be like this:

sol:= Optimization:-NLPSolve(
     f, {}, {(x1,x2,x3,x4)-> x1+x2+x3+x4-65}, 0..20, 0..15, 0..20, 0..15, minimize
);

The Units subpackages (Standard in this case) overload the arithmetic operators to behave like this. You can temporarily revert any operator to its default meaning by giving it the :- prefix and using it in prefix form. So, for this case, you could do

eq:= :-`+`(2.*Unit(m), d) = 5.*Unit(m);

     (the equation will display normally: the + will be infix.)

Now you can pass this equation to solve:

solve(eq, d);

e:= w^2/(w^4+2*w^2+1):
expand(numer(e)/w^2)/expand(denom(e)/w^2);

The first expand isn't really needed in this case.

Use 1. * Unit(m). It's not really meaningful to have an exact integer number of something with a unit.

The Hermite form of a matrix of polynomials is a good way to analyze the rank and how the rank depends on the values of the variables.

The following worksheet will not display on MaplePrimes, but please download it.

HermiteForm.mw

Here is the worksheet without output:



restart:

macro(LA= LinearAlgebra):

J1:= LA:-RandomMatrix((5,5)) + x*LA:-RandomMatrix((5,5), density= .15);

LA:-HermiteForm(J1);

So the rank of J depends on whether x is a root of the cubic polynomial in the 5,5 position. The rank can be 4 or 5.

 

Let's impose another level of rank deficiency:

Dr4:= LA:-DiagonalMatrix([1,1,1,0,0]):

J2:= expand~(J1.Dr4.LA:-RandomMatrix((5,5)));

expand~(LA:-HermiteForm(J2));

The rank is 3 independent of x.

J3:= expand~(J1.Dr4.J1^+);

H3:= expand~(LA:-HermiteForm(J3));

X:= solve(H3[3,3]=0);

LA:-ReducedRowEchelonForm(expand~(eval(J3, x= X)));

LA:-ReducedRowEchelonForm(expand~(eval(J3, x= 1)));

The rank is 3 independent of x, but the final reduced row-echelon form does depend on x.



Download HermiteForm_no_output.mw

By default, Maple carries 10 decimal digits through its floating-point computations. That's not enough to accurately evaluate f(31536000.).  You can change the number of digits of precision with

Digits:= 15;

for example.

From what I've been able to find out from some very cursory web searches, Mathematica uses "machine precision" by default. This corresponds to about 17 decimal digits (53 bits, to be precise). You can get Maple to use machine precision by setting

UseHardwareFloats:= true;

Maple has no practical limit on the number of Digits you can use. I believe that's true for Mathematica also.

From what I've been able to find out from very cursory web searches, the HP Prime uses and is limited to 12 digits, and the TI Nspire uses 10 by default and is limited to 14.

No matter what system you use or what digits setting you use, you will always be able to find some number for which this computation has significantly fewer accurate digits than you asked for. Indeed, it'll always be possible to find a large number x such that f(x) = 1.

*All information about Mathematica, TI, and HP is based on very quick ad hoc web searches. Please don't quote me on those.

 

To show that Maple knows the other half of the Fundamental Theorem of Calculus, I'd do

Diff(Int(f(t), t= a..x), x):  % = value(%);

First 218 219 220 221 222 223 224 Last Page 220 of 395