Carl Love

Carl Love

28020 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Here is a not-fully-polished implementation of Horner's method. It will find the non-multiple real roots of a polynomial with real coefficients. It will likely fail if it encounters a multiple root. Its output is the roots followed by the remaining unfactored polynomial.

Suggestions for improvement are most welcome.

Horner:= proc(P::polynom)
local
    p:= expand(P),
    x:= indets(p, And(name, Not(constant))),
    oldDigits:= Digits,
    roots:= table(),
    R, d
;    
    # Check that p is univariate.
     if nops(x) > 1 then
          error "Multiple names found: ", x
     else
          x:= x[]
     end if;
     
     Digits:= Digits+2+ilog10(Digits);

     # Get upper bound of abs val of roots by Rouche's theorem
     p:= evalf(collect(p,x));
     p:= expand(p/lcoeff(p));
     R:= 1+max(abs~([coeffs(p)]));
    
     for d to degree(p) do
         try
            R:= Student:-NumericalAnalysis:-Newton(
                 p, x= R, tolerance= 10^(-oldDigits-1), maxiterations= 2*Digits
            )
        catch:
            userinfo(1, Horner, StringTools:-FormatMessage(lastexception[2..-1]));
            break
        end try;
        p:= quo(p, x-R, x); # Reduce degree of p using the root
        roots[d]:= R        
    end do;
    
    evalf[oldDigits](sort(convert(roots, list)))[], evalf[oldDigits](p)
end proc:  

Example:

p:= `*`('randpoly(x, degree= 3)' $ 3);

infolevel[Horner]:= 1:
Horner(p);
Horner: maximum number of iterations (36) exceeded

fsolve(p);
 

 

Suppose that your Array is named A, and that its columns are indexed starting with 1 (the default). Then to plot the first vs. second columns

plot(A[.., [1,2]]); or plot(A[.., [2,1]]);

and to plot the first vs. third columns

plot(A[.., [1,3]]); or plot(A[.., [3,1]]);

exports(foo, static);

See ?exports .

There should be separate pages for each object and each object-method combination. Redundacy is acceptable in the "Description" sections, but try to vary the examples. Do not put things on the same page just because they are the same name overloaded.

This should be a Post, not a Question, because the answers are open-ended rather than definitive. Also, Posts tend to get more responses than Questions.

Yes, your commands are correct. But the answers returned by implicitdiff are more complicated than they need to be. After implicitdiff, use factor(simplify(%)).

See the examples at ?maximize or ?minimize .

The output of CrossProduct (or &x) is already a Vector, and does not need any conversion to be used with plots:-arrow.


with(LinearAlgebra):

v:= CrossProduct( < 1, 2, 3 > , < 4, 5, 6 > );

v := Vector(3, {(1) = -3, (2) = 6, (3) = -3})

plots:-arrow(v, axes= box);

 


Download arrow.mw

Using NonlinearFit is not appropriate in this case because there are three parameters in the Gompertz model and you only have three data points. I used solve instead.


restart:

Model:= t-> A*exp(B*exp(C*t));

proc (t) options operator, arrow; A*exp(B*exp(C*t)) end proc

Data:= [[0, 151326], [10, 179323], [20, 203202]]:

[{solve({seq(Model(d[1])=d[2], d= Data)}, {A,B,C})}[]]:

Sol:= evalf(%);

[{A = 212499.757668900, B = -.644068448786462, C = -.133346528689064}, {A = 288156.835274021, B = -.644068448786462, C = -0.305930714458130e-1}]

plot(
     [Data, seq(eval(Model(t), s), s= Sol)], t= -30..30,
     style= [point, line$2], color= [red, blue, green], symbolsize= 15
);

Obviously the blue curve does not satisfy the data (red diamonds). I don't know why solve includes it. We need to throw out Sol[1].

 

 

assign(Sol[2]);

plot(Model(t), t= -50..50);

 


Download Gompertz.mw

It works better if you just use inequalities:

restart:
assume(n >= 1);
is(n^2 <= 2^n);
                             false
assume(n >= 4);
is(n^2 <= 2^n);
                              true

Whenever is returns FAIL, you can retry the command with _EnvTry:= hard. In this case, it doesn't help.

You have used f in two different ways in your question---as the function name and as one of the coefficients. I changed the function name to g. Here's how to solve the problem:

restart:
g:= x-> d*x^2+e*x+f:
solve({g(0)=4, g(3)=7, g(-2)=18}, {d,e,f});

You need to use a Vector (or an Array) as temporary local storage for the list. Maple will not let you assign directly to a parameter.

InsertionSort:= proc(L::list)
local A:= < L >, j, key, i;
     for j from 2 to nops(L) do
          key:= A[j];
          for i from j by -1 to 2 while A[i-1] > key do
               A[i]:= A[i-1]
          end do;
          A[i]:= key
     end do;
     convert(A, list)
end proc;

Also note that the j loop completely surrounds the i loop, and that Vectors (and lists) are indexed starting at 1, not 0.

You usually cannot remove a removable singularity by using eval. You need to use limit instead. In your case, change

eval(diff(Psi_, epsilon), epsilon= 0);

to

limit(diff(Psi_, epsilon), epsilon= 0);

This will take a few more seconds to process, but definitely less than a minute.

Let T be a Vector of your time values, and let Y be a Vector of the corresponding population values (or whatever you are measuring if not population). Then, in its most basic form, the command would be

Statistics:-NonlinearFit(A*exp(B*exp(C*t)), T, Y, t);

See ?Statistics,NonlinearFit .

Usually, when you want to convert to binary, you want to access the individual bits/digits. The command Bits:-Split returns a list of the bits.

nums:= [10,3,90,6]:
Bits:-Split~(nums);
    [[0, 1, 0, 1], [1, 1], [0, 1, 0, 1, 1, 0, 1], [0, 1, 1]]

Note that each list of bits is returned in order from least-significant bit to most-significant bit. See ?Bits,Split .

The above pertains to converting nonnegative integers only.

To search for the global optimum in univariate problems with no general constraints but finite bounds, use method= branchandbound.

Note that the method is somewhat confused by the strict inequality at 10 in your piecewise. If you change it to x <= 10, then Maximize will return x = 10.

First 343 344 345 346 347 348 349 Last Page 345 of 395