Using method of eigenfunction expansion to solve a second order steady-state PDE. Does what I have look correct so far?

Question

I was following the process in Richard Haberman's Applied Partial Differential Equations, pgs. 388-389, to solve a PDE of the form $\nabla^2 u(x,y) = Q(x,y)$ which was subjected to homogeneous BCs using method of eigenfunction expansion. I tested my own example just to make sure and wanted to verify if what I had done to set up the ODE was correct. Here is the example I made:

$$u_{xx} + u_{yy} = 2\Big(x(x-1)+y(y-1)\Big); \ \ 0\le x,y \le 1 \\ BCs = \cases{u(0,y) = 0 & u(1,y) = 0 \\ u(x,0) = 0 & u(x,1) = 0}$$

We have homogeneous BCs going in both directions, and as I understand it this is to mean that we can use x-dependent or y-dependent eigenfunctions because both $\phi_n(y)$ and $\phi_n(x)$ satisfy the boundary conditions. The book decided to go with $\phi_n(x) = \sin(n\pi x/L)$ so the solution was written in the form: $$u_1 = \sum_{n=1}^\infty b_n(y) \sin(\frac{n\pi x}{L})$$

The next step was plugging this into the PDE yields the following ODE:

$$\sum_{n=1}^\infty \bigg[\frac{d^2 b_n}{dy^2}-\Big(\frac{n\pi}{L}\Big)^2 b_n\bigg]\sin(\frac{n\pi x}{L}) = Q(x,y)$$

We're going to utilize $\int_0^L \sin(n\pi x/L)\sin(m \pi x /L)dx = \cases{0 & m =/= n \\ L/2 & m=n}$ to eliminate the $\sum$ and $\sin(n \pi x /L)$ term and isolate $b_n(y)$. All we need to do is multiply both sides by $\sin(m \pi x /L)$ and integrate.

$$\int_0^L \Big(\sum_{n=1}^\infty \bigg[\frac{d^2 b_n}{dy^2}-\Big(\frac{n\pi}{L}\Big)^2 b_n\bigg]\sin(\frac{n\pi x}{L})\sin(\frac{m \pi x}{L})\Big)dx = \int_0^L \Big(Q(x,y)\sin(\frac{m \pi x}{L})\Big)dx \\ \sum_{n=1}^\infty \Big(\int_0^L \bigg[\frac{d^2 b_m(y)}{dy^2}-\Big(\frac{n\pi}{L}\Big)^2 b_m(y)\bigg]\sin(\frac{n\pi x}{L})\sin(\frac{m \pi x}{L})\Big)dx = \int_0^L \Big(Q(x,y)\sin(\frac{m \pi x}{L})\Big)dx \\ =0 + 0 + 0 + ... + \Big(\int_0^L \bigg[\frac{d^2 b_m(y)}{dy^2}-\Big(\frac{m\pi}{L}\Big)^2 b_m(y)\bigg]\sin^2(\frac{m\pi x}{L})\Big)dx+ 0 + 0 + 0 = \Big(Q(x,y)\sin(\frac{m \pi x}{L})\Big)dx \\ \boxed{\bigg[\frac{d^2 b_m(y)}{dy^2}-\Big(\frac{m\pi}{L}\Big)^2 b_m(y)\bigg] = \frac{2}{L}\int_0^L\Big(Q(x,y)\sin(\frac{m \pi x}{L})\Big)dx}$$

We'll shorthand $m\pi=k$ and then plug in the forcing term for Q(x,y), we get:

$$b(y)''-k^2b(y) = \frac{4}{L}\int_0^L\Big(x(x-1)+y(y-1)\Big)\sin(kx)dx \\ =4\bigg(\int_0^1\Big(x(x-1)\sin(kx)\Big)dx+\int_0^1\Big(y(y-1)\sin(kx)\Big)dx\bigg) \\ u = x(x-1), \ \ dv = \sin(kx) \\ du = x-1+x(1), \ \ v = -\frac{1}{k}\cos(kx) \\ =4\bigg(-\frac{1}{k}\cos(kx)x(x-1)|_0^1-(-\frac{1}{k}\int_0^1\Big((2x-1)\cos(kx)\Big)dx)-\frac{y(y-1)}{k}\cos(kx)|_0^1\bigg) \\ u = 2x-1, \ \ dv = \cos(kx) \\ du = 2, \ \ v = \frac{1}{k}\sin(kx) \\ =4\bigg(-(\frac{1}{k}\cos(k(1))1((1)-1)-\frac{1}{k}\cos(k(0))0((0)-1))+\frac{1}{k}\Big(\frac{1}{k}\sin(kx)(2x-1)|_0^1 - \frac{2}{k}\int_0^1\sin(kx)dx\Big)-\frac{y(y-1)}{k}(\cos(k(1))-\cos(k(0)))\bigg) \\ =4\bigg(0+\frac{1}{k}\Big(\frac{1}{k}\sin(kx)(2x-1)|_0^1 - \frac{2}{k}\int_0^1\sin(kx)dx\Big)-\frac{y(y-1)}{k}(\cos(k)-1)\bigg) \\ =4\bigg(\frac{1}{k}\Big(\frac{1}{k}\big(\sin(k(1))(2(1)-1)-\sin(k(0))(2(0)-1)\big) + \frac{2}{k}\big(\frac{1}{k}\cos(kx)|_0^1\big)\Big)-\frac{y(y-1)}{k}(\cos(k)-1)\bigg) \\ =4\bigg(\frac{1}{k}\Big(\frac{1}{k}\big(\sin(k)-0\big) + \frac{2}{k^2}\big((\cos(k(1))-\cos(k(0))\big)\Big)-\frac{y(y-1)}{k}(\cos(k)-1)\bigg) \\ =4\bigg(\frac{1}{k}\Big(\frac{1}{k}\sin(k) + \frac{2}{k^2}\big(\cos(k)-1\big)\Big)-\frac{y(y-1)}{k}(\cos(k)-1)\bigg) \\ \text{given k = } m\pi \text{ we know } \sin(k)=0\\ =4 \Big(y(y-1)\big(1-\cos(k)\big)-\frac{2}{k^2}\big(1-\cos(k)\big)\Big) \\ \text{and } \cos(k) = (-1)^m \\ \boxed{b(y)''-k^2b(y)=\frac{4}{m\pi}\Big(1-(-1)^m\Big)\bigg(y(y-1)-\frac{2}{(m\pi)^2}\bigg)}$$

And now we need to solve the second-order ODE. We start with the characteristic form to find the homogeneous solution for $b_m$ and then use some method to find the particular part.

$$r^2 - (m\pi)^2 = 0 \rightarrow r = \pm m\pi \\ b_g(y) = c_1 e^{-m\pi y}+c_2 e^{m\pi y} \\ \text{And then we need to find the particular part: } b_p = ?$$

Was the ODE (final boxed equation) set up correctly?

Answer 1

Too long for a comment.

I didn't check all your calculations --- you can use a computer algebra system to do that ---, but your method of solution is correct. I just want to call your attention to two other methods which, in my opinion, are simpler or faster:

1. The first one is similar to yours, but deals with $x$ and $y$ symmetrically: write the solution as $$ u(x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}A_{mn}\sin(m\pi x)\sin(n\pi y). \tag{1} $$ Plugging $(1)$ into the PDE, you obtain $$ u_{xx}+u_{yy}=-\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}(m^2+n^2)\pi^2A_{mn}\sin(m\pi x)\sin(n\pi y)=Q(x,y). \tag{2} $$ Now, use the orthogonality of the sines to obtain the coefficients $A_{mn}$: $$ A_{mn}=-\frac{4}{(m^2+n^2)\pi^2}\int_0^1\int_0^1Q(x,y)\sin(m\pi x)\sin(n\pi y)\,dxdy. \tag{3} $$
2. The second method only works in rare cases, but when it works it is very fast: guess the solution! My guess, in this case, is $u(x,y)=Ax(x-1)y(y-1)$. It satisfies all the boundary conditions, and $$ u_{xx}+u_{yy}=2Ay(y-1)+2Ax(x-1)=2A[x(x-1)+y(y-1)], \tag{4} $$ so it also satisfies the PDE if $A=1$.