Search code examples
maple

How to solve Laplace PDE on disk with piecewise boundary condition in Maple?

Maple 2017.3 on windows, is able to solve Laplace PDE on disk, when the boundary condition f(theta) is given as generic function (not defined).

But when I define f(theta) using piecewise, it does not solve the PDE for some reason. Any one knows why? Am I doing something wrong?

Here is MWE showing it can solve the PDE on disk of radius 1, with generic f(theta) for the boundary conditions

restart;
pde:= diff(u(r,theta),r$2)+1/r*diff(u(r,theta),r) + 1/r^2* diff(u(r,theta),theta$2)=0;
bc:= u(1,theta)=f(theta), 
     u(r,0)=u(r,2*Pi),
     D[2](u)(r,0)=D[2](u)(r,2*Pi);
sol:= pdsolve([pde,bc],u(r,theta)) assuming r<=1  and theta>=0 and theta<2*Pi;

Mathematica graphics

Now I changed f(theta) above, and gave it a specific function, which is piecewise. Now it does not solve it

restart;
f:=theta->piecewise(theta>=0 and theta<=Pi,20,
                    theta>Pi and theta<2*Pi,0);
pde:= diff(u(r,theta),r$2)+1/r*diff(u(r,theta),r) + 1/r^2* diff(u(r,theta),theta$2)=0;
bc:= u(1,theta)=f(theta), 
     u(r,0)=u(r,2*Pi),
     D[2](u)(r,0)=D[2](u)(r,2*Pi);
sol:= pdsolve([pde,bc],u(r,theta)) assuming r<=1 and theta>=0 and theta<2*Pi;

Mathematica graphics

Am I doing something wrong? Is there a work around?

For reference, Mathematica 11.2 can solve this:

ClearAll[u,r,θ]
pde=D[u[r,θ],{r,2}]+1/r D[u[r,θ],r]+1/r^2 D[u[r,θ],{θ,2}]==0
bc=Piecewise[{{20,0<θ<Pi},{0,Pi<θ<2Pi}}]
sol=DSolve[{pde,u[1,θ]==bc},u[r,θ],r,θ];
sol=sol/.K[1]->n

Mathematica graphics

Solution

  • So why not obtain the result from pdsolve for unassigned f, and then substitute the piecewise?

    restart;

pde:= diff(u(r,theta),r$2)+1/r*diff(u(r,theta),r) + 1/r^2* 
diff(u(r,theta),theta$2)=0:
bc:= u(1,theta)=f(theta), 
     u(r,0)=u(r,2*Pi),
     D[2](u)(r,0)=D[2](u)(r,2*Pi):

sol:= rhs(pdsolve([pde,bc],u(r,theta)))
      assuming r<=1 and theta>=0 and theta<2*Pi:

lprint(sol);
   (1/2)*(2*(Sum(r^n*((Int(f(theta)*sin(n*theta), theta = 0 .. 
   2*Pi))*sin(n*theta)+(Int(f(theta)*cos(n*theta), theta = 0 .. 
   2*Pi))*cos(n*theta))/Pi, n = 1 .. infinity))*Pi+Int(f(theta), theta = 0 .. 2*Pi))/Pi

F:=theta->piecewise(theta>=0 and theta<=Pi,20,
                    theta>Pi and theta<2*Pi,0):

    Below we'll use subs(f=(F),sol) for the pde solution with the piecewise F substituted for the unassigned name f.

    Let's compute sol at a specific point.

    sol_spec:=simplify(combine(eval(subs(f=(F),sol),[r=1/2,theta=Pi/3]))):

value(sol_spec): # symbolic summation and symbolic integration

lprint(%);
   (1/2)*((80/3)*Pi+40*arctan((1/5)*3^(1/2)))/Pi

evalf(%);
                      15.45628948

evalf(sol_spec); # numeric summation and numeric integration
                      15.51328895

    We can resolve the integrals in sol and obtain the summation form (as claimed for Mma).

    T:=subs(f=F,subsindets(sol,specfunc(Int),IntegrationTools:-Split,Pi)):
ans:=simplify(eval(T,Int=int)):

lprint(%);
   10+Sum(20*sin(n*theta)*(r^n-(-r)^n)/(n*Pi), n = 1 .. infinity)

    Evaluate that at the same specific point as before.

    ans_spec:=(eval(ans,[r=1/2,theta=Pi/3])):

value(ans_spec): # symbolic summation

lprint(%);
   40/3+20*arctan((1/5)*3^(1/2))/Pi

evalf(%);
                      15.45628948

evalf(ans_spec); # numeric summation
                      15.51328895

    We can even obtain a closed form for the summation under some conditions on r and theta.

    huh:=evalc(subsindets(ans,specfunc(Sum),
                      u->sum(op(u),formal)))
     assuming r>0, r<1, theta>0, theta<Pi:

lprint(huh);
   10+20*arctan(sin(theta)*r/(1-cos(theta)*r))/Pi+20*arctan(sin(theta)*r/(1+cos(theta)*r))/Pi

eval(huh,[r=1/2,theta=Pi/3]):

lprint(%);
   40/3+20*arctan((1/5)*3^(1/2))/Pi

evalf(%);
                      15.45628948

    Note that there is a bug (already reported) in the combining of pairs of symbolic arctan call, which prevents me from simplifying or combining them in huh above. This doesn't affect the above results.