Carl Love

Carl Love

28050 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

The output data standarderrors is only available for a linear fit. You can easily convert your problem to a linear one. 

X1:= ln~(X); Y1:= ln~(Y -~ 1);

Now Y1 = ln(a) + b*X1.

Fit(a1+b*n, X1, Y1, n, output= [...]);

Finally a = exp(a1).

Did you have in mind a procedure to plot the circles and the eigenvalues?

restart:

Gershgorin:= proc(M::Matrix(square,complexcons))
uses LA= LinearAlgebra, P= plots;
local A:= evalf(M), n:= LA:-RowDimension(A), i, j;
     P:-display(
          [
               seq(
                    plottools:-circle(
                         [Re,Im](A[i,i]),
                         add(abs(A[i,j]), j= 1..n) - abs(A[i,i]),
                         color= COLOR(HSV, i/n, .5, .5)
                    ),
               i= 1..n
               ),
               P:-pointplot(
                    [Re,Im]~(LA:-Eigenvalues(A)),
                    color= red
               )
          ],
          scaling= constrained,
          symbol= diagonalcross, symbolsize= 16,
          _rest          
     )
end proc:

Example:

macro(LA= LinearAlgebra):

A:= LA:-RandomMatrix(5)+I*LA:-RandomMatrix(5);

A := Matrix(5, 5, {(1, 1) = 29-89*I, (1, 2) = -14+96*I, (1, 3) = 44+34*I, (1, 4) = 11-32*I, (1, 5) = 27-81*I, (2, 1) = 35+95*I, (2, 2) = 37+69*I, (2, 3) = 92-55*I, (2, 4) = 61-9*I, (2, 5) = 58+11*I, (3, 1) = -70+77*I, (3, 2) = -97+72*I, (3, 3) = 73+54*I, (3, 4) = 28+69*I, (3, 5) = 2-76*I, (4, 1) = -43-84*I, (4, 2) = -92+42*I, (4, 3) = -39+79*I, (4, 4) = -48+31*I, (4, 5) = 54+82*I, (5, 1) = -23-63*I, (5, 2) = 73+55*I, (5, 3) = 62-99*I, (5, 4) = -63-66*I, (5, 5) = 47-29*I})

Gershgorin(A);

A:= LA:-RandomMatrix(5) + I*LA:-RandomMatrix(5)+
    7*LA:-DiagonalMatrix(
         LA:-RandomVector(5)+I*LA:-RandomVector(5)
    );

A := Matrix(5, 5, {(1, 1) = 42+32*I, (1, 2) = 82-63*I, (1, 3) = 12+12*I, (1, 4) = 22+21*I, (1, 5) = 60-82*I, (2, 1) = -32+91*I, (2, 2) = 107-453*I, (2, 3) = -62+45*I, (2, 4) = 14+90*I, (2, 5) = -95-70*I, (3, 1) = -1-I, (3, 2) = 42+30*I, (3, 3) = -670-196*I, (3, 4) = 16+80*I, (3, 5) = -20+41*I, (4, 1) = 52+63*I, (4, 2) = 18+10*I, (4, 3) = -68+60*I, (4, 4) = -299-121*I, (4, 5) = -25+91*I, (5, 1) = -13-23*I, (5, 2) = -59+22*I, (5, 3) = -67-35*I, (5, 4) = 99+88*I, (5, 5) = -215-517*I})

Gershgorin(A);

 

 

Download Gershgorin.mw

 

 

You may be able to achieve what you want using the command ?plots,textplot . If you could provide a picture of what you want, I'll try to make it in Maple.

Perhaps you are referring to the tick marks on the axes? It is quite easy to change the position of them. See ?plot,tickmarks .

Don't use capital I as a variable. It is reserved for the imaginary unit. I usually use J instead when I want to use I as a variable.

Don't assign r or d.

Then the algsubs command will work.

What does "y = -0,5+7" mean? Do you mean y= -0.5*x+7? In that case,

plot([2*x-3, -0.5*x+7]);

Here is a not-fully-polished implementation of Muller's method. It will likely fail on multiple roots, and there is no checking for division by zero.

restart:

Muller:= proc(
     f::appliable, inits::[complexcons,complexcons,complexcons],
     {tolerance:= 10.^(1-Digits)}, {maxsteps:= 99}
)
local
     d1, d2, d, b, D, p, k,
     oldDigits:= Digits,
     x0:= inits[1], x1:= inits[2], x2:= inits[3],
     h1:= x1-x0, h2, h,
     fx0, fx1, fx2
 ;
     Digits:= Digits+2+ilog10(Digits);
     fx0:= evalf(f(x0));  fx1:= evalf(f(x1));
     d1:= (fx1-fx0)/h1;
     for k to maxsteps do
          h2:= x2-x1;
          fx2:= evalf(f(x2));
          d2:= (fx2-fx1)/h2;
          d:= (d2-d1)/(h2+h1);
          b:= d2+h2*d;
          D:= evalf(sqrt(b^2-4*fx2*d));
          h:= 2*fx2/(b+`if`(abs(b-D) < abs(b+d), D, -D));
          p:= x2-h;
          userinfo(2, Muller, p);
          if abs(h) < tolerance then
               userinfo(1, Muller, sprintf("Used %d steps", k));  
               return evalf[oldDigits](p)  
          end if;
          x0:= x1;  x1:= x2;  x2:= p;
          fx0:= f(x0);  fx1:= fx2;  fx2:= f(p);
          h1:= h2;  d1:= d2;
     end do;
     FAIL
end proc:
          
infolevel[Muller]:= 2:
Muller(x-> exp(x)+1, [1, 0, -1]);
Muller: -1.08197670686932642-1.58493557680537191*I
Muller: -1.27350123007187934-2.85054027345107362*I
Muller: -.381407414149884648-3.74751157829847389*I
Muller: .197272946056780758-3.14572035519676954*I
Muller: -0.16781896044693894e-1-3.15707732044571707*I
Muller: 0.9078666240418e-6-3.14209510958315899*I
Muller: 0.260490388101188484e-6-3.14159295025712521*I
Muller: -0.80434247308e-13-3.14159265359054070*I
Muller: -0.2472981852860e-18-3.14159265358979324*I
Muller: -0.247298185265482997e-18-3.14159265358979324*I
Muller: Used 10 steps
                         

Check:

evalf(Pi);

                                             3.14159265358979

Compare with fsolve:

fsolve(exp(x)+1, x=-I, complex);
                            

@MOSO1401 Sure, I can help you to make an animation. What do you want to vary? I guess that the above plot represents t=0 and you want to show how the wave evolves over time. Here is the most basic technique:

N:= 20:
for k from 1 to N do
    t:= k*Pi/N;
    frame[k]:= plot3d(sin(x+t)*cos(y+t), x= -Pi..Pi, y= -Pi..Pi)
end do:
plots:-display([seq(frame[k], k= 1..N)], insequence);

F:= eval(f(t), dsnx):
G:= eval(g(t), dsnx):

Now F and G are real-valued procedures of real arguments. For example,

F(1);
                  0.32274425109407956

If you want, you can also extract the data matrices directly from the plot. Let me know if that is what you meant.

plots:-surfdata(IS);

Everything that you want will happen if you give the command

with(RealDomain);

rem(6*x^2+5*x+1, 2*x+1, x);
                               0

The second polynomial must be a factor of the first:

factor(6*x^2+5*x+1);
                      (2 x + 1) (3 x + 1)

Maple 17's integration is more sophisticated with respect to using assumptions.

int((1+(a*x)^t)^(-1/t), x= 0..z) assuming t>0, t<1;

It would be very difficult to get your expression into exactly the form that you want. But here are some easy tricks that get it fairly close.


Expression entered without ln.

f:= product((p*beta[1]*(t[i]/theta[1])^(beta[1]-1)*exp(-(t[i]/theta[1])^beta[1])/theta[1])^Y[i]*beta[2]*(t[i]/theta[2])^beta[2]/(t[i]*exp((t[i]/theta[2])^beta[2])*((1-p)*beta[2]*(t[i]/theta[2])^(beta[2]-1)*exp(-(t[i]/theta[2])^beta[2])/theta[2])^Y[i])-(p*beta[1]*(t[i]/theta[1])^(beta[1]-1)*exp(-(t[i]/theta[1])^beta[1])/theta[1])^Y[i]*beta[2]*(t[i]/theta[2])^beta[2]*p/(t[i]*exp((t[i]/theta[2])^beta[2])*((1-p)*beta[2]*(t[i]/theta[2])^(beta[2]-1)*exp(-(t[i]/theta[2])^beta[2])/theta[2])^Y[i]), i = 1 .. n);

product((p*beta[1]*(t[i]/theta[1])^(beta[1]-1)*exp(-(t[i]/theta[1])^beta[1])/theta[1])^Y[i]*beta[2]*(t[i]/theta[2])^beta[2]/(t[i]*exp((t[i]/theta[2])^beta[2])*((1-p)*beta[2]*(t[i]/theta[2])^(beta[2]-1)*exp(-(t[i]/theta[2])^beta[2])/theta[2])^Y[i])-(p*beta[1]*(t[i]/theta[1])^(beta[1]-1)*exp(-(t[i]/theta[1])^beta[1])/theta[1])^Y[i]*beta[2]*(t[i]/theta[2])^beta[2]*p/(t[i]*exp((t[i]/theta[2])^beta[2])*((1-p)*beta[2]*(t[i]/theta[2])^(beta[2]-1)*exp(-(t[i]/theta[2])^beta[2])/theta[2])^Y[i]), i = 1 .. n)

I use Int instead of sum below because the Int can be distributed with IntegrationTools:-Expand. I'll convert it back to sum later.

applyop(ln, 1, subs(product= Int, f));

Int(ln((p*beta[1]*(t[i]/theta[1])^(beta[1]-1)*exp(-(t[i]/theta[1])^beta[1])/theta[1])^Y[i]*beta[2]*(t[i]/theta[2])^beta[2]/(t[i]*exp((t[i]/theta[2])^beta[2])*((1-p)*beta[2]*(t[i]/theta[2])^(beta[2]-1)*exp(-(t[i]/theta[2])^beta[2])/theta[2])^Y[i])-(p*beta[1]*(t[i]/theta[1])^(beta[1]-1)*exp(-(t[i]/theta[1])^beta[1])/theta[1])^Y[i]*beta[2]*(t[i]/theta[2])^beta[2]*p/(t[i]*exp((t[i]/theta[2])^beta[2])*((1-p)*beta[2]*(t[i]/theta[2])^(beta[2]-1)*exp(-(t[i]/theta[2])^beta[2])/theta[2])^Y[i])), i = 1 .. n)

simplify(%) assuming positive;

Int(-t[i]^beta[1]*theta[1]^(-beta[1])*Y[i]+ln(t[i])*Y[i]*beta[1]-ln(t[i])*Y[i]*beta[2]-ln(theta[1])*Y[i]*beta[1]+beta[2]*Y[i]*ln(theta[2])+Y[i]*t[i]^beta[2]*theta[2]^(-beta[2])+beta[2]*ln(t[i])-beta[2]*ln(theta[2])-t[i]^beta[2]*theta[2]^(-beta[2])+ln(beta[1])*Y[i]-Y[i]*ln(beta[2])-ln(t[i])+ln((1-p)^(-Y[i])*p^Y[i]-(1-p)^(-Y[i])*p^(1+Y[i]))+ln(beta[2]), i = 1 .. n)

simplify(%, symbolic);

Int(-t[i]^beta[1]*theta[1]^(-beta[1])*Y[i]+ln(t[i])*Y[i]*beta[1]-ln(t[i])*Y[i]*beta[2]-ln(theta[1])*Y[i]*beta[1]+beta[2]*Y[i]*ln(theta[2])+Y[i]*t[i]^beta[2]*theta[2]^(-beta[2])+beta[2]*ln(t[i])-beta[2]*ln(theta[2])-t[i]^beta[2]*theta[2]^(-beta[2])+ln(beta[1])*Y[i]-Y[i]*ln(beta[2])-ln(t[i])-Y[i]*ln(1-p)+ln(p^Y[i]-p^(1+Y[i]))+ln(beta[2]), i = 1 .. n)

IntegrationTools:-Expand(%);

-(Int(t[i]^beta[1]*Y[i], i = 1 .. n))/theta[1]^beta[1]+beta[1]*(Int(ln(t[i])*Y[i], i = 1 .. n))-beta[2]*(Int(ln(t[i])*Y[i], i = 1 .. n))-beta[1]*ln(theta[1])*(Int(Y[i], i = 1 .. n))+beta[2]*ln(theta[2])*(Int(Y[i], i = 1 .. n))+(Int(Y[i]*t[i]^beta[2], i = 1 .. n))/theta[2]^beta[2]+beta[2]*(Int(ln(t[i]), i = 1 .. n))-beta[2]*ln(theta[2])*(Int(1, i = 1 .. n))-(Int(t[i]^beta[2], i = 1 .. n))/theta[2]^beta[2]+ln(beta[1])*(Int(Y[i], i = 1 .. n))-ln(beta[2])*(Int(Y[i], i = 1 .. n))-(Int(ln(t[i]), i = 1 .. n))-ln(1-p)*(Int(Y[i], i = 1 .. n))+Int(ln(p^Y[i]-p^(1+Y[i])), i = 1 .. n)+Int(ln(beta[2]), i = 1 .. n)

eval(%, Int= sum);

-n*beta[2]*ln(theta[2])+n*ln(beta[2])+beta[2]*ln(theta[2])*(sum(Y[i], i = 1 .. n))-beta[1]*ln(theta[1])*(sum(Y[i], i = 1 .. n))+beta[1]*(sum(ln(t[i])*Y[i], i = 1 .. n))-beta[2]*(sum(ln(t[i])*Y[i], i = 1 .. n))+(sum(Y[i]*t[i]^beta[2], i = 1 .. n))/theta[2]^beta[2]-(sum(t[i]^beta[1]*Y[i], i = 1 .. n))/theta[1]^beta[1]+ln(beta[1])*(sum(Y[i], i = 1 .. n))+beta[2]*(sum(ln(t[i]), i = 1 .. n))-(sum(t[i]^beta[2], i = 1 .. n))/theta[2]^beta[2]-ln(beta[2])*(sum(Y[i], i = 1 .. n))-ln(1-p)*(sum(Y[i], i = 1 .. n))+sum(ln(p^Y[i]-p^(1+Y[i])), i = 1 .. n)-(sum(ln(t[i]), i = 1 .. n))


Download ln_of_product.mw

See ?readdata . It's like ImportMatrix, but it lets you specify a separate datatype for each column.

To get the complex roots of a polynomial p(z), use fsolve(p(z), complex). The Maple command for the argument is simply argument. The modulus is the same as the absolute value, which is abs. The maximum value in a list can be found with max.

Below, I also included a plot of the roots of the two polynomials in the complex plane.

R1:= [fsolve(z^5-2*z^4-2*z^3-2, complex)];

[-.856360086695669-.527505295612770*I, -.856360086695669+.527505295612770*I, .476701432159350-.699464811611423*I, .476701432159350+.699464811611423*I, 2.75931730907264]

argument~(R1);

[-2.58950203832279, 2.58950203832279, -.972578893785131, .972578893785131, 0.]

max(%);

2.58950203832279

R2:= [fsolve(z^5-5*z^3-2*z-2, complex)];

[-2.28375944914536, -.569085775991655, .252021816767795-.769192677734267*I, .252021816767795+.769192677734267*I, 2.34880159160143]

abs~(R2);

[2.28375944914536, .569085775991655, .809427187341117, .809427187341117, 2.34880159160143]

max(%);

2.34880159160143

 

plot(
     [[Re,Im]~(R1), [Re,Im]~(R2)],
     style=point, symbolsize= 20, symbol= [cross,diamond],
     axes= box, scaling= constrained, view= [-3..3,-1..1]
);

 

 

Download arg_and_modulus.mw

First 341 342 343 344 345 346 347 Last Page 343 of 395