Robert Israel

6577 Reputation

21 Badges

18 years, 212 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are answers submitted by Robert Israel

That printlevel := -1;  threw me for a while.  Don't do that!

Unfortunately, point doesn't seem to like a indexed names.  You could try something like

>  centpoint:= [seq](point(CentPoint||i,cent[i]),i=1..20);

Actually, let's remove the Re, which prevents the expression from being analytic.  I'll also use exact rational arithmetic rather than floats.

> G := 0.11111e-5*((-50000.*(-1.0182*R^18+8.9640*10^22-6.5799*10^21*R^3+4.0243*10^19*R^6+1.6412*10^18*R^9+13416.*sqrt(-1.*R^6*(3.6362*10^28*R^12+1.8565*10^10*R^21-7.9260*10^30*R^9-6.9380*10^18*R^15+5.9390*10^32*R^6+35821.*R^24+7.4270*10^16*R^18+1.8038*10^35-1.7654*10^34*R^3)))^(2/3)-1.0491*10^17*R^6-50601.*R^12+4.9005*10^18*R^3+2.0120*10^5*R^6*(-1.0182*R^18+8.9640*10^22-6.5799*10^21*R^3+4.0243*10^19*R^6+1.6412*10^18*R^9+13416.*sqrt(-1.*R^6*(3.6362*10^28*R^12+1.8565*10^10*R^21-7.9260*10^30*R^9-6.9380*10^18*R^15+5.9390*10^32*R^6+35821.*R^24+7.4270*10^16*R^18+1.8038*10^35-1.7654*10^34*R^3)))^(1/3)-1.0950*10^11*R^3*(-1.0182*R^18+8.9640*10^22-6.5799*10^21*R^3+4.0243*10^19*R^6+1.6412*10^18*R^9+13416.*sqrt(-1.*R^6*(3.6362*10^28*R^12+1.8565*10^10*R^21-7.9260*10^30*R^9-6.9380*10^18*R^15+5.9390*10^32*R^6+35821.*R^24+7.4270*10^16*R^18+1.8038*10^35-1.7654*10^34*R^3)))^(1/3)+4.4754*10^12*(-1.0182*R^18+8.9640*10^22-6.5799*10^21*R^3+4.0243*10^19*R^6+1.6412*10^18*R^9+13416.*sqrt(-1.*R^6*(3.6362*10^28*R^12+1.8565*10^10*R^21-7.9260*10^30*R^9-6.9380*10^18*R^15+5.9390*10^32*R^6+35821.*R^24+7.4270*10^16*R^18+1.8038*10^35-1.7654*10^34*R^3)))^(1/3)+(86605.*I)*(-1.0182*R^18+8.9640*10^22-6.5799*10^21*R^3+4.0243*10^19*R^6+1.6412*10^18*R^9+13416.*sqrt(-1.*R^6*(3.6362*10^28*R^12+1.8565*10^10*R^21-7.9260*10^30*R^9-6.9380*10^18*R^15+5.9390*10^32*R^6+35821.*R^24+7.4270*10^16*R^18+1.8038*10^35-1.7654*10^34*R^3)))^(2/3)-(1.8171*10^17*I)*R^6-(87648.*I)*R^12+(8.4882*10^18*I)*R^3-1.0015*10^20-1.7347*10^20*I)/(R^6*(-1.0182*R^18+8.9640*10^22-6.5799*10^21*R^3+4.0243*10^19*R^6+1.6412*10^18*R^9+13416.*sqrt(-1.*R^6*(3.6362*10^28*R^12+1.8565*10^10*R^21-7.9260*10^30*R^9-6.9380*10^18*R^15+5.9390*10^32*R^6+35821.*R^24+7.4270*10^16*R^18+1.8038*10^35-1.7654*10^34*R^3)))^(1/3)));
 

> GR:= convert(G, rational);
   A:=asympt(GR, R, 12);

The result is a rather complicated expression involving the csgn function.  It should be possible to get rid of that in the limit R -> infinity.

 

It should be possible to get rid of that csgn function as R -> infinity.

> cs:= op(indets(A, specfunc(anything, csgn)));

cs := csgn(1/5000*I*(-5091+448200000000000000000000000/R^18-32899500000000000000000000/R^15+201215000000000000000000/R^12+8206000000000000000000/R^9+67080000*(-36362000000000000000000000000/R^12-18565000000/R^3+7926000000000000000000000000000/R^15+6938000000000000000/R^9-593900000000000000000000000000000/R^18-35821-74270000000000000/R^6-180380000000000000000000000000000000/R^24+17654000000000000000000000000000000/R^21)^(1/2)/R^3))

> A:= eval(A, cs = csgn(limit(op(cs),R=infinity)));

A := 1/900009*(-50000*exp(2/3*ln(5091/5000)+2/3*I*Pi)-50601-87648*I+201200*exp(1/3*ln(5091/5000)+1/3*I*Pi)+86605*I*exp(2/3*ln(5091/5000)+2/3*I*Pi))*exp(-1/3*ln(5091/5000)-1/3*I*Pi)+(22360000/4581945819*(-50000*exp(2/3*ln(5091/5000)+2/3*I*Pi)-50601-87648*I+201200*exp(1/3*ln(5091/5000)+1/3*I*Pi)+86605*I*exp(2/3*ln(5091/5000)+2/3*I*Pi))*exp(-1/3*ln(5091/5000)-1/3*I*Pi)*exp(1/2*ln(35821)+1/2*I*Pi)+1/900009*(2236000000000/5091*exp(2/3*ln(5091/5000)+2/3*I*Pi)*exp(1/2*ln(35821)+1/2*I*Pi)-4498832000000/5091*exp(1/3*ln(5091/5000)+1/3*I*Pi)*exp(1/2*ln(35821)+1/2*I*Pi)-3872975600000/5091*I*exp(2/3*ln(5091/5000)+2/3*I*Pi)*exp(1/2*ln(35821)+1/2*I*Pi)-109500000000*exp(1/3*ln(5091/5000)+1/3*I*Pi))*exp(-1/3*ln(5091/5000)-1/3*I*Pi))/R^3+O(1/(R^6))

Now back to floating point:

> evalf(A);

-.4915365171e-5-.1108660606e-5*I+(-121662.2429+5.456694018*I)/R^3+O(1/(R^6))

 

 

My guess is that you have somehow set up an infinite recursion where some entry of your Matrix is defined in terms of itself. 
But it's just a guess.  Why don't you show us your Matrix?

The problem is that you're using e.  This just another variable to Maple, not the constant 2.71....  You need to use the exponential function exp.

Thus:

> int(exp(-x^2/2), x = -infinity .. infinity);

 

Can you give an example of such a procedure?  In what form are you saving it?

ProbabilityFunction seems not to be very clever: it can handle a random variable with a named (discrete) distribution, but not something more complicated.  However, you can use Probability.

> with(Statistics): Y := RandomVariable(Geometric(1/p)) + 1;
> Probability(Y = i) assuming i::posint;

1/p*(1-1/p)^(-1+i)

Maple can't quite do the integral by itself. It needs a bit of help.

> S:=expand(ln(x^(alpha[1]-alpha[2])*(1-x)^(beta[1]-beta[2])*Beta(alpha[2],beta[2])
/Beta(beta[1],alpha[1]))) assuming x>0,x<1, alpha[1]>0,alpha[2]>0,beta[1]>0,beta[2]>0;


 

 

S := (alpha[1]-alpha[2])*ln(x)+(beta[1]-beta[2])*ln(1-x)+ln(Beta(alpha[2],beta[2])/Beta(alpha[1],beta[1]))

> J:= map(t -> Int(x^(alpha[1]-1)*(1-x)^(beta[1]-1)*t,x=0..1)/Beta(alpha[1],beta[1]),S);
J := 1/Beta(alpha[1],beta[1])*Int(x^(alpha[1]-1)*(1-x)^(beta[1]-1)*(alpha[1]-alpha[2])*ln(x),x = 0 .. 1)+1/Beta(alpha[1],beta[1])*Int(x^(alpha[1]-1)*(1-x)^(beta[1]-1)*(beta[1]-beta[2])*ln(1-x),x = 0 .. 1)+1/Beta(alpha[1],beta[1])*Int(x^(alpha[1]-1)*(1-x)^(beta[1]-1)*ln(Beta(alpha[2],beta[2])/Beta(alpha[1],beta[1])),x = 0 .. 1)
> Jv:= map(value,J) assuming alpha[1]>0,beta[1]>0;

Jv := (alpha[1]-alpha[2])*(Beta(alpha[1],beta[1])*Psi(alpha[1])-Beta(alpha[1],beta[1])*Psi(alpha[1]+beta[1]))/Beta(alpha[1],beta[1])+1/Beta(alpha[1],beta[1])*int(x^(alpha[1]-1)*(1-x)^(beta[1]-1)*(beta[1]-beta[2])*ln(1-x),x = 0 .. 1)+ln(Beta(alpha[2],beta[2])/Beta(alpha[1],beta[1]))

> J1:= op(indets(Jv,specfunc(anything,int)));

J1 := int(x^(alpha[1]-1)*(1-x)^(beta[1]-1)*(beta[1]-beta[2])*ln(1-x),x = 0 .. 1)

> IntegrationTools[Change](subs(int=Int,J1),x=1-t);

(beta[1]-beta[2])*Int((1-t)^(alpha[1]-1)*t^(beta[1]-1)*ln(t),t = 0 .. 1)

> V1:=value(%) assuming alpha[1]>0,beta[1]>0;

V1 := (beta[1]-beta[2])*(Beta(alpha[1],beta[1])*Psi(beta[1])-Beta(alpha[1],beta[1])*Psi(alpha[1]+beta[1]))

> Jv:= subs(J1=V1,Jv);

Jv := (alpha[1]-alpha[2])*(Beta(alpha[1],beta[1])*Psi(alpha[1])-Beta(alpha[1],beta[1])*Psi(alpha[1]+beta[1]))/Beta(alpha[1],beta[1])+(beta[1]-beta[2])*(Beta(alpha[1],beta[1])*Psi(beta[1])-Beta(alpha[1],beta[1])*Psi(alpha[1]+beta[1]))/Beta(alpha[1],beta[1])+ln(Beta(alpha[2],beta[2])/Beta(alpha[1],beta[1]))

> collect(Jv,Psi);

(alpha[1]-alpha[2])*Psi(alpha[1])+(beta[1]-beta[2])*Psi(beta[1])+(-alpha[1]+alpha[2]-beta[1]+beta[2])*Psi(alpha[1]+beta[1])+ln(Beta(alpha[2],beta[2])/Beta(alpha[1],beta[1]))

 

The problem seems to be that you are using Zeta, which for Maple is not just a Greek letter but the Riemann zeta function.  If you change Zeta to zeta, it should work.

> tf1:= subs(Zeta = zeta, tf);
  evalc(Re(tf1));

-M[tot]/R[eq]*(-omega^2*M[tot]/R[eq]-2*omega^2*zeta*omega[n]*C[p]+omega[n]^2*M[tot]/R[eq])/((-omega^2*M[tot]/R[eq]-2*omega^2*zeta*omega[n]*C[p]+omega[n]^2*M[tot]/R[eq])^2+(-M[tot]*C[p]*omega^3+omega*omega[n]^2*M[tot]*C[p]+2*omega*zeta*omega[n]*M[tot]/R[eq]-omega*Theta^2)^2)-M[tot]*C[P]*omega*(-M[tot]*C[p]*omega^3+omega*omega[n]^2*M[tot]*C[p]+2*omega*zeta*omega[n]*M[tot]/R[eq]-omega*Theta^2)/((-omega^2*M[tot]/R[eq]-2*omega^2*zeta*omega[n]*C[p]+omega[n]^2*M[tot]/R[eq])^2+(-M[tot]*C[p]*omega^3+omega*omega[n]^2*M[tot]*C[p]+2*omega*zeta*omega[n]*M[tot]/R[eq]-omega*Theta^2)^2)

> evalc(Im(tf1));

-M[tot]*C[P]*omega*(-omega^2*M[tot]/R[eq]-2*omega^2*zeta*omega[n]*C[p]+omega[n]^2*M[tot]/R[eq])/((-omega^2*M[tot]/R[eq]-2*omega^2*zeta*omega[n]*C[p]+omega[n]^2*M[tot]/R[eq])^2+(-M[tot]*C[p]*omega^3+omega*omega[n]^2*M[tot]*C[p]+2*omega*zeta*omega[n]*M[tot]/R[eq]-omega*Theta^2)^2)+M[tot]/R[eq]*(-M[tot]*C[p]*omega^3+omega*omega[n]^2*M[tot]*C[p]+2*omega*zeta*omega[n]*M[tot]/R[eq]-omega*Theta^2)/((-omega^2*M[tot]/R[eq]-2*omega^2*zeta*omega[n]*C[p]+omega[n]^2*M[tot]/R[eq])^2+(-M[tot]*C[p]*omega^3+omega*omega[n]^2*M[tot]*C[p]+2*omega*zeta*omega[n]*M[tot]/R[eq]-omega*Theta^2)^2)

 

 

You must have done something else besides what you wrote, because I get

> o:= map(e -> e^z mod n, g);

o := [1, 2, 3]

Either g was not [1, 128, 130], or z was not 23, or n was not 187, when you executed the command that gave you [0,0,0].

It's hard to tell for sure, because you're still posting pictures instead of text, but I suspect that what you have is e.g.

    r^2*cos^2*theta

instead of

   r^2*cos(theta)^2

cos and sin are functions, and their arguments need to be inside parentheses, despite the common typesetting convention that allows omitting the parentheses.

So you want

> f:= (r, theta) -> r/(2*sqrt(r^2*cos(theta)^2 + r^2*sin(theta)^2));
   evalf(int(int(f(r,theta),r=1/sin(2*theta)*cos(theta)..2*sin(theta)),theta=Pi/8..Pi/2));

.5201568037

In Maple 12 I only get the plot going up to value 1, not 4500.  This behaviour
seems to have changed between Maple 11 and Maple 12.  But note that you specified range = 0 .. 0.1, so you shouldn't be surprised if strange things happen outside that range.  Try it with range = 0 .. 1 (the help file doesn't make it clear, but the "range" should be for the variable that is not "time").

 

Your X2 and X1 are the same.  This won't produce very good results.  But anyway:

> with(Statistics):
   f:= a+b*x[1]+c*x[2];
   X:= Matrix([X1,X2])^%T;
   Fit(f, X,Y,[x[1],x[2]]);

Warning, model is not of full rank
3.90909090909090740+.590909090909089496*x[1]+.590909090909092605*x[2]

 

 

Since your objective will be discontinuous, it's probably not a good idea to formulate it this way.  You might try
three separate problems, one for 2*x[1]^2+3*x[2]+4*x[3]+2*x[4] <= 10, one for 10 <= 2*x[1]^2+3*x[2]+4*x[3]+2*x[4]<= 15,
and one for 2*x[1]^2+3*x[2]+4*x[3]+2*x[4] >= 15.

Something like this:

> PDE:=diff(p(a,t),t)+diff(p(a,t),a)=-1;
> IBC:={p(a,0)=sin(a),p(0,t)=cos(t)};
> pds:=pdsolve(PDE,IBC,numeric,time=t,range=0..10); 
   Q:=subs(pds:-value(t=2,output=listprocedure),p(a,t));
   evalf(Int(Q,0..10));

                     -16.23669314

(a rather inaccurate result, as the correct answer is sin(2)-17-cos(8) = -15.94520254)

sprintf("%.4f", 0.1234);

 

First 76 77 78 79 80 81 82 Last Page 78 of 138