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

Some general principles to make numerical analysis routines run faster are

  1. set Digits to 15 (in general, but not always),
  2. create all Matrices and Vectors with option datatype= float[8],
  3. use evalhf as much as possible.

See ?evalhf and ?LinearAlgebra,General,EffNumLA .

You should post some code that seems to run particularly slowly, and I'll show you how to improve its efficiency.

I don't think that it makes sense to maximize one's chances (chances of winning presumably). I think that the thing to maximize is the expected value of participating.

The optimal strategy is going to be the same for everybody, and there's a sort of prisoner's dilemma involved here. The optimal straegy is going to be something akin to this: The employees who want to participate hold a preliminary drawing to select three of their members. Then those three, and only those three, each buy one ticket. Then those 3 have a 100% chance of winning. If that is not feasible, then each employee should use a computer to draw a random integer from 1 to 1600. Anyone who draws a 1, 2, or 3 should buy one ticket. Then each of those has close to a 100% chance of winning.

Is this what you had in mind?

plots:-fieldplot([piecewise(x+y > 0, x), piecewise(x+y > 0, y)], x= -10..10, y= -10..10);

 

Let x represent the overall expression, let n represent the numerator, and let d represent the denominator. We have the following system of three equations:

solve({x=1+n/d, n=3+x/n, d=2+d/x}, {x,n,d});

Note that iquo and irem can be done together with a single call to either. Here's how I'd write your procedure:

euclid:= proc(m::posint, n::posint)
local q, r1:= m, r2:= n, r3;
     while r2 <> 0 do
          q:= iquo(r1,r2,'r3');  #or r3:= irem(r1,r2,'q');
          printf("%d = %d*%d + %d\n", r1, q, r2, r3);
          (r1,r2):= (r2,r3)
     end do;
     r1
end;

To get f(6) and f(15), just enter the procedure into Maple, then enter f(6) and f(15).

All four can be entered directly into dsolve. The first two require the addition of option implicit


Problem 1:

restart:

dsolve({(x+epsilon*y(x))*diff(y(x),x)+y(x)=0, y(1)=1}, implicit);

-(1+(1/2)*epsilon)/y(x)+x+(1/2)*epsilon*y(x) = 0

solve(%, y(x));

(-x+(epsilon^2+x^2+2*epsilon)^(1/2))/epsilon, -(x+(epsilon^2+x^2+2*epsilon)^(1/2))/epsilon

Problem 2:

restart:

dsolve({(x^n+epsilon*y(x))*diff(y(x),x)+n*x^(n-1)*y(x)=m*x^(m-1), y(1)=b}, implicit);

y(x)*x^n-x^m+(1/2)*epsilon*y(x)^2-b+1-(1/2)*epsilon*b^2 = 0

solve(%, y(x));

(-x+(b^2*epsilon^2+2*x^m*epsilon+2*b*epsilon+x^2-2*epsilon)^(1/2))/epsilon, -(x+(b^2*epsilon^2+2*x^m*epsilon+2*b*epsilon+x^2-2*epsilon)^(1/2))/epsilon

Problem 3:

restart:

Sol:= dsolve({diff(u(t),t$2)+u(t)+epsilon*u(t)^3=0, u(0)=A, D(u)(0)=0});

u(t) = RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))*((-2*epsilon-4)/((epsilon*RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))^2-epsilon-2)*(epsilon+2)))^(1/2)*JacobiSN(((1/2)*(2*epsilon+4)^(1/2)*t+InverseJacobiSN(A/(RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))*(-2/(epsilon*RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))^2-epsilon-2))^(1/2)), RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))*(-(epsilon+2)*epsilon)^(1/2)/(epsilon+2))/(-2/(epsilon*RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))^2-epsilon-2))^(1/2))*((-2*epsilon-4)/((epsilon*RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))^2-epsilon-2)*(epsilon+2)))^(1/2), RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))*(-(epsilon+2)*epsilon)^(1/2)/(epsilon+2))

indets(%, specfunc(anything, RootOf));

{RootOf(_Z*2^(1/2)*((A^2*_Z^2*epsilon-A^2*epsilon-2*A^2+2*_Z^2)/_Z^2)^(1/2)*(-(2*A^2*_Z^2*epsilon^2-2*A^2*epsilon^2-4*A^2*epsilon-4*epsilon-8)/(epsilon+2))^(1/2))}

alias(R=%[]);

R

Sol;

u(t) = R*((-2*epsilon-4)/((epsilon*R^2-epsilon-2)*(epsilon+2)))^(1/2)*JacobiSN(((1/2)*(2*epsilon+4)^(1/2)*t+InverseJacobiSN(A/(R*(-2/(epsilon*R^2-epsilon-2))^(1/2)), R*(-(epsilon+2)*epsilon)^(1/2)/(epsilon+2))/(-2/(epsilon*R^2-epsilon-2))^(1/2))*((-2*epsilon-4)/((epsilon*R^2-epsilon-2)*(epsilon+2)))^(1/2), R*(-(epsilon+2)*epsilon)^(1/2)/(epsilon+2))

Problem 4:

restart:

dsolve({diff(y(x),x$2)+2*diff(y(x),x)^2+3*diff(y(x),x) = -sin(x), y(0)=1, D(y)(0)=0});

y(x) = -(3/4)*x+(1/2)*ln((MathieuC(-9, -4, -(1/4)*Pi+(1/2)*x)*MathieuSPrime(-9, -4, (1/4)*Pi)+3*MathieuC(-9, -4, -(1/4)*Pi+(1/2)*x)*MathieuS(-9, -4, (1/4)*Pi)+3*MathieuS(-9, -4, -(1/4)*Pi+(1/2)*x)*MathieuC(-9, -4, (1/4)*Pi)+MathieuS(-9, -4, -(1/4)*Pi+(1/2)*x)*MathieuCPrime(-9, -4, (1/4)*Pi))/(MathieuSPrime(-9, -4, -(1/4)*Pi+(1/2)*x)*MathieuC(-9, -4, -(1/4)*Pi+(1/2)*x)-MathieuCPrime(-9, -4, -(1/4)*Pi+(1/2)*x)*MathieuS(-9, -4, -(1/4)*Pi+(1/2)*x)))+1

 


Download 4_dsolves.mw

Let X:= [x1, ..., xm] and Y:= [[y11, ..., y1n], ..., [ym1, ..., ymn]]. Then

plot([[[X[i],Y[i][j]] $ i= 1..m] $ j= 1..n], style= point);

This version uses color to show which Y-list the point came from.

Let me know how that goes.

Here are some hints: Each geometric object asked for can be created with a single command from ?geom3d with no intermediate computation required. Just follow the syntax shown in the help files. The commands required are, in order: with, ?geom3d,point , point, point, ?geom3d,line , ?geom3d,plane , ?geom3d,intersection , ?geom3d,sphere (with the the centername option),  line, ?geom3d,FindAngle , ?geom3d,coordinates , and ?geom3d,intersection . And you'll need evalf to get the decimal approximation for the first answer.

You just need to set the environment variable ?UseHardwareFloats :

UseHardwareFloats:= false;

Let me know how that goes.

(You probably know this, but I'll say it anyway: Maple's floating-point system is base 10, not base 2. So don't expect some of the standard examples of floating-point errors to work exactly the same way as is shown in the textbooks.)

You need to put forget(x) before the x:= proc....

Since Maple allows a newline character anywhere that it allows whitespace, continuation characters are very rarely needed. Your example would work just as well without the backslash. Even long strings can be split "simply "
"like "
"this."

There is no logical concept of a "line" of code in Maple. Statements can, and often do, go on for tens of lines.

That being said, I have no information about the situation in Maple IDE.

Instead of type(e, negative) you need to say is(e < 0). The type of an expression is an inherent property of the expression, and it is not affected by assumptions.

You should not use and in the assumptions; just separating them with the commas is logically equivalent to and.

Putting it together, your statement should be

is(e < 0) assuming  
     d1 > 0, d2 > 0, d3 > 0, Kw > 0, Ks > 0,
     0 < pfw1, pfw1 < 1, 0 < pfw2, pfw2 < 1, 0 < pfw3, pfw3 < 1,
     0 < pfs1, pfs1 < 1, 0 < pfs2, pfs2 < 1, 0 < pfs3, pfs3 < 1,
     d1-pfw1*Kw-pfs1*Ks > 1, d2-pfw2*Kw-pfs2*Ks > 1, d3-pfw3*Kw-pfs3*Ks > 1,
     pfw1 <> pfw2, pfw2 <> pfw3, pfs1 <> pfs2, pfs2 <> pfs3
;

Maple's response is false.

Addressing your question 2: There is no algorithm for deriving the conditions under which a complicated, non-polynomial, multivariate expression is negative (or positive, or zero).

Let's say that your expression with the h for the difference quotient is called dq. Then you need to say

limit(dq, h= 0);

Here's your problem:

f:= x-> x^3+3*cos(x):
dq:= (f(x+h)-f(x))/h:
limit(dq, h=0);

 

Change plot3d to plots:-implicitplot3d. Then adjust the ranges of x, y, and z to get a better plot.

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