Robert Israel

6447 Reputation

21 Badges

15 years, 145 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are answers submitted by Robert Israel

pdsolve can sometimes find "symbolic" solutions for certain types of PDE's with boundary conditions, but it's rather hit-or-miss, and mostly miss.  In any case, boundary conditions involving a limit are certainly much too subtle for pdsolve.  You might try pdsolve without the boundary conditions at x=0 and x -> infinity to get ingredients for a separation-of-variables solution.

> pdsolve({diff(V(x,y),x$2)+diff(V(x,y),y$2)=0, V(x,0)=0,V(x,Pi)=0}, HINT=X(x)*Y(y));

Since you want the solution to go to 0 as x -> infinity, you'd want to take only exponentials of negative coefficients times x: in general

where the coefficients a[n] are found by expanding V(0,y) in a sine Fourier series.

In your case V(0,y) = sin(y), so you just need one term in this series:

10 billion points?  Yes, that's rather a lot.  Somewhat better if you use combinations rather than permutations (assuming "649" is a reference to the 6/49 lottery, the order of selected numbers is immaterial).  So there are only binomial(49,6) = 13983816 choices.  Now

combinat:-choose(49,6)

produces an "object too large" error.  So we'll have to be cleverer.

> with(ImageTools): 
M:= 1000; N:= 1000; # image width and height
    R:= Create(M, N,1);
    xmax:= binomial(49,6); ymax:= 444546474849;
    x:= 1;
    for a from 1 to 44 do
       for b from a+1 to 45 do
          for c from b+1 to 46 do
             for d from c+1 to 47 do
for e from d+1 to 48 do
for f from e+1 to 49 do
y:= parse(cat(a,b,c,d,e,f));
i:= round(1 + (x-1)/(xmax-1)*(M-1));
j:= round(1 + (y-1)/(ymax-1)*(N-1));
R[i,j]:= 1;
x:=x+1
od od od od od od:
Write("649.jpg",R);

 

Yes, it's a bug, probably related to difficulties solve has with trig functions.  The simplest thing to do is to express the CDF for Y "by hand" using the CDF for X:

> CDFY:= unapply( CDF(X,arctan(t)),t);

and then differentiate that to get the PDF.

> PDFY:= D(CDFY);

The functionality you're looking for is provided by casesplit in the PDEtools package.

For example,

> PDEtools[casesplit]([ diff(a(x)^2*b(x), x) = a(x)^2, diff(a(x)*b(x),x) = b(x)], [a(x), b(x)]);

 

It would have been better if you could have uploaded a Maple worksheet, so I could have tried this with your equations, rather than a .doc file.  I really don't feel like trying to input your complicated formulas into Maple by hand.

implicitplot3d is a very resource-intensive way to plot a sphere.  You might use sphere in the plottools package.  So e.g. something like:

> with(plots): with(plottools):
   f:= t -> display(sphere([0,0,0],0.5, grid=[15,10]),
                    seq(sphere([i*cos(t/i),i*sin(t/i),0],0.3, grid=[10,6]),i=1..9));
   animate(f,[t],t=0..20, scaling=constrained);

 

 

I see three possible issues.

  • 10*tan(x)^alpha contains two variables, x and alpha.  I assume you want this to be a function of x, in which case you'll want to either assign alpha a numerical value or evaluate this at a numerical value for alpha (e.g. using subs or eval).
  • tan(x) has negative values in some intervals.  If alpha is not an even integer, tan(x)^alpha will have complex values.  These will not show up on your plot.
  • tan(x) has a singularity at odd multiples of Pi/2.  So tan(x) will have very large values when x is near Pi/2.  If this is within the interval of x values for your plot, you'll want to restrict the y interval to see anything reasonable.

However, you seem to be interested in the interval 0 to Pi/4, for which the second and third issues are not relevant. 

This has been discussed here a number of times.  My favourite technique is

> f:= x^2-1: g:= -x-1;
plottools:-transform(unapply([x,y+g],x,y))(plot(f-g,x=-1.5 .. 1.5,filled=true));

Which equations are you trying to solve?  You have, for example, eqn1 which is

and eqn1b which is

You might be able to have one of these, but obviously not both.

Something like this?

> with(plots):

   animate(arrow,[[0,0],[0,cos(t)]],t=0..4*Pi,view=[-1..1,-1..1]);

Your TrapRule procedure is expecting a function, not an expression, as its first argument.  You

could try TrapRule(x -> exp(x), 4, -3, 1)

or just TrapRule(exp, 4, -3, 1);

I assume you want p p instead of theta in that solve command (since theta doesn't appear in the equations).  I see you already used factor.  You're not going to get any more factorization than that unless you specify a field extension.  In this case I guessed that there might be some factorization over the field generated by I, and I was right, although it doesn't make your expression look simpler.

> factor(rsoln, I);

Certainly there are no factors such as x1^2+x3^2+y1^2+y3^2.

Maybe a bit less obvious:

F:= n -> Matrix(n, n, charfcn[n+1] @ `+`);

Do you  mean something like this (for a positive integer n)?

Matrix(n,n, (i,j) -> piecewise(i-j mod n = 1, -1/2, i-j mod n = n -1, 1/2, 0));

or

Matrix(n,n,shape=Circulant[[0,1/2,0$n-3,-1/2]]);

@edgar : By default, Maple looks for a factorization over the field generated by the coefficients.  If you want a factorization over some algebraic number field, you can provide that as a second argument to factor, e.g.

> factor(x^2 - 2, sqrt(2));

For a univariate polynomial, you can also ask for a complete factorization over the (floating-point) real or complex field, e.g.

> factor(x^5 - 3*x + 2, real);

> factor(x^5 - 3*x + 2, complex);

 

>  expr:= diff(f(x,y),x);
    PDEtools:-dchange({x=r*cos(theta),y=r*sin(theta),f(x,y)=F(r,theta)},diff(f(x,y),x));

 

1 2 3 4 5 6 7 Last Page 2 of 138