stokes_paradox.mws

>    restart;

We express (Laplacian)^2 of psi in polar coordinates.

>    Lpsi := diff(psi(r,theta),r$2)+1/r*diff(psi(r,theta),r) + 1/r^2*diff(psi(r,theta),theta$2);
L2psi := diff(Lpsi,r$2)+1/r*diff(Lpsi,r) + 1/r^2*diff(Lpsi,theta$2):
collect(L2psi,r);

We compute (Laplacian)^2 of the nth term of the Fourier series (in theta) of psi.

>    nterm := psi[n](r)*exp(I*n*theta);
Lnterm := diff(nterm,r$2)+1/r*diff(nterm,r) + 1/r^2*diff(nterm,theta$2);
L2nterm := collect(diff(Lnterm,r$2)+1/r*diff(Lnterm,r) + 1/r^2*diff(Lnterm,theta$2),r):
collect(simplify(L2nterm*exp(-I*n*theta)),r)*exp(I*n*theta);

nterm := psi[n](r)*exp(n*theta*I)

Lnterm := diff(psi[n](r),`$`(r,2))*exp(n*theta*I)+1/r*diff(psi[n](r),r)*exp(n*theta*I)-1/r^2*psi[n](r)*n^2*exp(n*theta*I)

(diff(psi[n](r),`$`(r,4))+2*diff(psi[n](r),`$`(r,3))/r+(-diff(psi[n](r),`$`(r,2))-2*diff(psi[n](r),`$`(r,2))*n^2)/r^2+(diff(psi[n](r),r)+2*diff(psi[n](r),r)*n^2)/r^3+(-4*psi[n](r)*n^2+psi[n](r)*n^4)/r^...

We solve the ODE that results from setting (Laplacian)^2 of the nth term of the Fourier series equal to zero.

>    ode1 := diff(psi[n](r),r$4) + 2/r*diff(psi[n](r),r$3) - (2*n^2+1)/r^2*diff(psi[n](r),r$2)
         + (2*n^2+1)/r^3*diff(psi[n](r),r) + (n^4 - 4*n^2)/r^4*psi[n](r) = 0;
dsolve(ode1);

ode1 := diff(psi[n](r),`$`(r,4))+2*diff(psi[n](r),`$`(r,3))/r-(1+2*n^2)/r^2*diff(psi[n](r),`$`(r,2))+(1+2*n^2)/r^3*diff(psi[n](r),r)+(n^4-4*n^2)/r^4*psi[n](r) = 0

psi[n](r) = _C1*r^(n+2)+_C2*r^(-n)+_C3*r^n+_C4*r^(-n+2)

But the solution has a different form when n = 0, 1 and -1.

>    for j from 0 to 1 do
sprintf("for n = %d:",j);
ode1 := diff(psi[j](r),r$4) + 2/r*diff(psi[j](r),r$3) - (2*j^2+1)/r^2*diff(psi[j](r),r$2)
         + (2*j^2+1)/r^3*diff(psi[j](r),r) + (j^4 - 4*j^2)/r^4*psi[j](r) = 0;
dsolve(ode1);
end do;

ode1 := diff(psi[0](r),`$`(r,4))+2/r*diff(psi[0](r),`$`(r,3))-1/r^2*diff(psi[0](r),`$`(r,2))+1/r^3*diff(psi[0](r),r) = 0

psi[0](r) = _C1*r^2+_C2*r^2*ln(r)+_C3+_C4*ln(r)

ode1 := diff(psi[1](r),`$`(r,4))+2/r*diff(psi[1](r),`$`(r,3))-3/r^2*diff(psi[1](r),`$`(r,2))+3/r^3*diff(psi[1](r),r)-3/r^4*psi[1](r) = 0

psi[1](r) = _C1/r+_C2*r^3+_C3*r+_C4*r*ln(r)

So these are the general solutions:

>    psi[n] := r -> _C1*r^(-n+2)+_C2*r^(n+2)+_C3*r^(-n)+_C4*r^n;
psi[0] := r -> _C1*r^2+_C2*r^2*ln(r)+_C3+_C4*ln(r);
psi[1] := r -> _C1*r^3+_C2/r+_C3*r+_C4*r*ln(r);

psi[n] := proc (r) options operator, arrow; _C1*r^(-n+2)+_C2*r^(n+2)+_C3*r^(-n)+_C4*r^n end proc

psi[0] := proc (r) options operator, arrow; _C1*r^2+_C2*r^2*ln(r)+_C3+_C4*ln(r) end proc

psi[1] := proc (r) options operator, arrow; _C1*r^3+_C2/r+_C3*r+_C4*r*ln(r) end proc

>   

Maple TM is a registered trademark of Waterloo Maple Inc.
Math rendered by WebEQ