Carl Love

Carl Love

28015 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

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

Just make the existing code into a procedure with parameters mu, rho, and g. Like this:

restart: Digits:=20:

ThreePlots:= proc(mu, rho, g)
local
     C_0:= 0.4, d:= 0.0002,
     rho_v:= 1000.0, m:= evalf((4/3)*Pi*rho_v*(d/2)^3),   

     F_D:=proc(v::list)                 
     local
          vx:= v[1], vy:= v[2],
          k:= -evalf((3*Pi*mu*d+(Pi/8)*C_0*rho*(d^2)*sqrt(vx^2+vy^2))/m)
     ;      
          [k*vx,k*vy]
     end proc,                                                      

     Eskritt:=proc(vo::list,h)  
     local F:= F_D(vo);
          [vo[1]+h*F[1], vo[2]+h*(F[2]-g)]
     end proc,

     v_0:= 1.0, h:= 0.002,
 
     OnePlot:= proc(phi, N, C)
     local
          v0:=evalf([v_0*cos(phi), v_0*sin(phi)]),
          i, j, v, r
     ;
          v[0]:=v0;
          for i to N do  v[i]:= Eskritt(v[i-1],h)  end do:
          r[0]:=[0,0]:
          for i to N do  r[i]:= r[i-1]+h*v[i-1]  end do:
          plot([[r[j] $j=0..N]], color= C)
      end proc,

      p1:= OnePlot(Pi/4, 65, green),
      p2:= OnePlot(Pi/5, 55, red),
      p3:= OnePlot(Pi/6, 50, blue)
 ;     
     plots:-display(
          [p1,p2,p3], labels=["",""], scaling= constrained
     )
end proc:

ThreePlots(1.8*10^(-5), 1.2, 9.81);

Then call it again with the new values of mu, rho, and g.

You need to specify a value of alpha in the same places as you specified M. I used alpha=1, but the resulting plots were not within your view. You should eliminate the view option until the plot is otherwise as you want.

b:= A[..,3]:
C:= A[..,[1,2,4..11]]:
x:= LinearAlgebra:-LinearSolve(C,b):
x[8];

You need to define the function f like this:

f:= (t,v)-> 0.0207*v^2 - 893.58;

and not like this:

f(t,v):= ....

Perhaps you are confused by the ?piecewise notation in which the integral answer is presented. You can convert it back to abs form, and you can also convert it to ?Heaviside form.

 

restart:

int(abs(x-2), x);

piecewise(x <= 2, -(1/2)*x^2+2*x, 2 < x, (1/2)*x^2-2*x+4)

lprint(%);

piecewise(x <= 2, -(1/2)*x^2+2*x, 2 < x, (1/2)*x^2-2*x+4)

convert(%, abs);

(1/2)*abs(x-2)*x-abs(x-2)+2

convert(%, Heaviside);

Heaviside(x-2)*x^2-4*Heaviside(x-2)*x-(1/2)*x^2+2*x+4*Heaviside(x-2)

 

 

Download convert_abs.mw

The Interp command that I showed for one of your previous questions works just as well in the two-variable case. So here's your problem (iii):

restart:
U:= [0,1,2]:
V:= [x^2+7, x^3+2*x+3, x^3+5]:
Interp(U,V,y) mod 11;

To see the algorithm for doing this, enter

showstat(`mod/Interp/Interp`);

You just need to pay attention to lines 11-26.

G:= GF(7, 3, x^3+x+1):
G:-`/`(G:-ConvertIn(x^2));

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