$\nabla^{4}f=1$ for $f$ defined in a square from $x\in[-1,1]$ and $y\in[-1,1]$
The boundary conditions are:
$f=0,f_{xx}=0$ on $x=\pm 1$
$f=0,f_{yy}=\mp 1$ on $y=\pm 1$
I intend to solve this problem numerically but would like to find a series solution to check the accuracy of my scheme. Does one know how to find a series solution?
The differential operator is $$ \nabla^4 = D_x^4 + 2D_y^2D_x^2 + D_y^4. $$
Step 1: Find eigenvectors of $D_x^4$ satisfying the boundary conditions in the $x$-direction
We choose the $x$-direction because the boundary conditions are homogeneous.
To find eigenfunctions with respect to $x$, suppose $u$ is a function of $x$ satisfying $\nabla^4 u = \lambda^4 u$ for some $\lambda$. (I could have used $\lambda$ instead of $\lambda^4$. This choice just makes future expressions prettier.)
It follows that $$ D_x^4 u - \lambda^4 u = 0. $$
If $\lambda = 0$, the solution to this ODE is $$ u = c_3 x^3 + c_2 x^2 + c_1 x + c_0. $$ Using the boundary conditions in the $x$-direction $u(-1) = u(1) = u''(-1) = u''(1) = 0$, it is easy to verify that $c_0 = c_1 = c_2 = c_3 = 0$. That means there is no eigenfunction with respect to $\lambda = 0$. (An eigenvector must be non-zero by definition.)
Next, consider the case $\lambda \ne 0$. The solution of $D_x^4 u - \lambda^4 u = 0$ is $$ u = c_1 \cos \lambda x + c_2 \sin \lambda x + c_3 \cosh \lambda x + c_4 \sinh \lambda x. $$ Compute $u''$: $$ u'' = -c_1 \lambda^2 \cos \lambda x - c_2 \lambda^2 \sin \lambda x + c_3 \lambda^2 \cosh \lambda x + c_4 \lambda^2 \sinh \lambda x. $$ Use the boundary conditions in the $x$ direction: $u(-1) = u(1) = u''(-1) = u''(1) = 0$. This gives a system of equations \begin{align} c_1 \cos\lambda - c_2 \sin\lambda + c_3\cosh\lambda - c_4\sinh\lambda & = 0 \tag 1\\ c_1 \cos\lambda + c_2 \sin\lambda + c_3\cosh\lambda + c_4\sinh\lambda & = 0 \tag 2\\ -c_1 \cos\lambda + c_2 \sin\lambda + c_3\cosh\lambda - c_4\sinh\lambda & = 0 \tag 3\\ -c_1 \cos\lambda - c_2 \sin\lambda + c_3\cosh\lambda + c_4\sinh\lambda & = 0. \tag 4 \end{align} $(1) + (2) + (3) + (4)$ gives $c_3\cosh\lambda = 0$. Since $\cosh\lambda > 0$ for all $\lambda$, $c_3 = 0$.
$(1) - (2) + (3) - (4)$ gives $c_4\sinh \lambda = 0$. Since $\sinh\lambda = 0$ only when $\lambda = 0$, we get $c_4 = 0$.
$u$ becomes $u = c_1\cos\lambda x + c_2\sin\lambda x$, and the four equations reduce to two: \begin{align} c_1 \cos \lambda - c_2 \sin \lambda & = 0 \tag 5\\ c_1 \cos \lambda + c_2 \sin \lambda & = 0. \tag 6 \end{align} Considering $c_1$ and $c_2$ as unknowns, this system of linear equations would have a non-trivial solution if and only if $$ \det\begin{pmatrix} \cos\lambda & -\sin\lambda \\ \cos\lambda & \sin\lambda \end{pmatrix} = 0, $$ which is true if and only if $\cos\lambda \sin\lambda = 0$. This means $\cos\lambda = 0$ or $\sin\lambda = 0$.
Let us consider the case $\cos\lambda = 0$ first. This amounts to $\lambda$ being of the form $n\pi - \frac{\pi}2$, $n \in \mathbb N$. $(5)$ (or $(6)$) gives $c_2 = 0$. $u$ becomes $u = c_1 \cos \lambda x$. This means $\cos \lambda x$ is an eigenfunction corresponding to eigenvalue $\lambda^4$. (Eigenvectors can be scaled by a non-zero scalar, so we simply set $c_2 = 1$.)
The other case when $\sin\lambda = 0$ is similar. The possible values of $\lambda$ are $n\pi$, $n \in \mathbb N$. $(5)$ (or $(6)$) gives $c_1 = 0$, so $u = c_2\sin\lambda x$. The eigenfunction is $\sin\lambda x$.
For ease of notation, define $\lambda_n = n\pi - \frac{\pi}2$, $u_n = \cos \lambda_n x$, $\mu_n = n\pi$ and $v_n = \sin \mu_n x$. It follows that $(\lambda_n^4, u_n)$ and $(\mu_n^4, v_n)$ are eigenvalue-eigenvector pairs of $D_x^4$ with respect to the boundary conditions in the $x$-direction. (In this setting, it is common to hear people say that $\lambda_n$ and $\mu_n$ are eigenvalues instead of $\lambda_n^4$ and $\mu_n^4$.)
We have discovered ALL eigenvalue-eigenfunction pairs of $D_x^4$ satisfying the homogeneous boundary conditions in the $x$-direction.
Note the following relations between eigenfunctions and related differential operators: \begin{align} D_x^4 u_n & = \lambda_n^4 u_n \tag 7\\ D_x^4 v_n & = \mu_n^4 v_n \tag 8\\ D_x^2 u_n & = -\lambda_n^2 u_n \tag 9\\ D_x^2 v_n & = -\mu_n^2 v_n \tag{10} \end{align} These will be needed later.
Step 2: Assume an expression of $f$ as a combination of eigenfunctions and compute $\nabla^4 f$
For a fixed $y$, the function $x \mapsto f(x, y)$ is a function in $x$ satisfying the boundary conditions in the $x$-direction. We assume that it can be expressed as a combination of eigenfunctions. More precisely, assume
$$ f(x, y) = \sum_{n=1}^n \left( a_n(y)u_n(x) + b_n(y)v_n(x) \right). \tag{11} $$ I'll drop the parentheses and assume $f$, $a_n$, $b_n$, $u_n$ and $v_n$ are all functions of appropriate variables. With the help of $(7), (8), (9), (10)$, we can compute \begin{align} \nabla^4 f & = (D_x^4 + 2D_y^2D_x^2 + D_y^4)f \\ & = \sum_{n=1}^\infty \left( \lambda_n^4 a_n u_n + \mu_n^4 b_n v_n - 2\lambda_n^2 (D_y^2 a_n) u_n - 2\mu_n^2 (D_y^2 b_n) v_n + (D_y^4 a_n) u_n + (D_y^4 b_n) v_n \right) \\ & = \sum_{n=1}^\infty \left( u_n \cdot \left(\lambda_n^4 - 2\lambda_n^2 D_y^2 + D_y^4\right)a_n + v_n \cdot \left(\lambda_n^4 - 2\lambda_n^2 D_y^2 + D_y^4\right)b_n \right). \tag{12} \end{align}
Step 3: Express the non-homogeneous forcing term as a combination of eigenfunctions
We're solving $\nabla^4 f = g$, and we have an expression for $f$ as above. We would like to have an expression for $g$ that looks similar to that of $f$. I'll describe the general method for an arbitrary $g$ first. The case $g(x, y) = 1$ is computed at the end of this step.
Since the expansion is in $x$-direction, it suffices to discuss the case where $g$ is only a function of $x$ first. If $g$ can be written as $$ g = \sum_{n=1}^\infty (c_n u_n + d_n v_n) $$ where $c_n$ and $d_n$ are constants, then by orthogonality of Fourier basis (recall that $u_n = \cos\left(\left(n - \frac 12\right)\pi x\right)$ and $v_n = \sin\left(n\pi x\right)$), we get \begin{align} c_n & = \frac{\langle g, u_n \rangle}{\langle u_n, u_n\rangle} = \int_{-1}^1 g(x)u_n(x) dx \tag{13}\\ d_n & = \frac{\langle g, v_n \rangle}{\langle v_n, v_n\rangle} = \int_{-1}^1 g(x)v_n(x) dx \tag{14} \end{align} Now the general case where $g$ is a function of $x$ and $y$ can be handled just by doing the same thing for each $y$. More precisely, \begin{align} g(x, y) & = \sum_{n=1}^\infty\left(c_n(y)u_n(x) + d_n(y)v_n(x)\right) \tag{15}\\ c_n(y) & = \frac{\langle g, u_n \rangle}{\langle u_n, u_n\rangle} = \int_{-1}^1 g(x, y)u_n(x) dx \tag{16}\\ d_n(y) & = \frac{\langle g, v_n \rangle}{\langle v_n, v_n\rangle} = \int_{-1}^1 g(x, y)v_n(x) dx \tag{17} \end{align}
The special case $g = 1$ gives \begin{align} c_n & = \int_{-1}^1 \cos\left(\left(n - \frac 12\right)\pi x\right) dx = \frac{1}{\left(n - \frac 12\right)\pi} = \frac{1}{\lambda_n} \tag{18}\\ d_n & = \int_{-1}^1 \sin\left(\left(n - \frac 12\right)\pi x\right) dx = 0 \tag{19} \end{align} (Please check my math.)
Step 4: Express non-homogeneous boundary conditions in terms of eigenfunctions
Similar to $g$, the boundary conditions in $y$ direction should be written in terms of eigenfunctions of $D_x^4$. (Each boundary condition is a function of $x$.)
Let $f_{-1}(x) = f(x, -1)$, $f_1(x) = f(x, 1)$, $F_{-1}(x) = f_{yy}(x, -1)$ and $F_1(x) = f_{yy}(x, 1)$. The given problem has $f_{-1} = f_1 = 0$, $F_{-1} = 1$ and $F_1 = -1$. $f_{-1}$ and $f_1$ are obviously zero, so no computation is needed. Since $F_{-1} = g$, we get the same coefficients of $F_{-1}$ as those of $g$. $F_1 = -g$, so the coefficients are just the opposite of those of $g$.
Apply these conditions to $f$ to get \begin{align} 0 = f_{-1} & = \sum_{n=1}^\infty \left(u_n \cdot a_n(-1) + v_n \cdot b_n(-1) \right) \\ 0 = f_{1} & = \sum_{n=1}^\infty \left(u_n \cdot a_n(1) + v_n \cdot b_n(1) \right) \\ D_y^2 f & = \sum_{n=1}^\infty \left(u_n \cdot D_y^2 a_n + v_n \cdot D_y^2 b_n\right) \\ g = F_{-1} & = \sum_{n=1}^\infty \left(u_n \cdot a''_n(-1) + v_n \cdot b''_n(-1)\right) \\ -g = F_{1} & = \sum_{n=1}^\infty \left(u_n \cdot a''_n(1) + v_n \cdot b''_n(1)\right) \\ \end{align} Comparing coefficients (which is justified by orthogonality of Fourier basis) gives \begin{align} a_n(-1) & = b_n(-1) = 0 \tag{I1}\\ a_n(1) & = b_n(1) = 0 \tag{I2}\\ a''_n(-1) & = c_n = \frac{1}{\lambda_n}\tag{I3}\\ b''_n(-1) & = d_n = 0 \tag{I4}\\ a''_n(1) & = -c_n = -\frac{1}{\lambda_n}\tag{I5}\\ b''_n(1) & = -d_n = 0 \tag{I6}\\ \end{align}
Step 5: Compare coefficients of both sides of the PDE at each $y$
Now we solve $\nabla^4 f = g$ by comparing coefficients of in front of $u_n$ and $v_n$ at each $y$. This method of comparing coefficients is justified by orthogonality of Fourier basis. From $(12)$ and $(15)$, \begin{align} (\lambda_n^4 - 2\lambda_n^2 D_y^2 + D_y^4)a_n & = c_n \tag{20}\\ (\mu_n^4 - 2\mu_n^2 D_y^2 + D_y^4)b_n & = d_n. \tag{21} \end{align} These two equations are ODEs in one variable, $y$. The differential operator in both equations are of the same form:
$$ (D_y^4 - 2\nu^2 D_y^2 + \nu^4)p = q $$ where $p$ and $q$ are functions of $y$. The characteristic polynomial, with $\alpha$ as an indeterminate, is $\alpha^4 - 2\nu^2\alpha^2 + \nu^4 = (\alpha^2 - \nu^2)^2 = (\alpha - \nu)^2(\alpha + \nu)^2$. Therefore, the null-space of the differential operator is spanned by $\{\cosh\nu y, y\cosh\nu y, \sinh\nu y, y\sinh\nu y\}$. In the case of a general $q$, $p$ can be obtained by variation of parameters. The solution will be of the form $p(y) = p_1\cosh\nu y + p_2y\cosh\nu y + p_3\sinh\nu y + p_4y\sinh\nu y + \int K(y, \tilde y) q(\tilde y) d\tilde y$ for some function $K$.
When $\nu$, $p$ and $q$ are replaced by terms in $(20)$ or $(21)$, we get \begin{align} a_n = & \alpha_n\cosh\lambda_n y + \beta_n y\cosh\lambda_n y + \gamma_n\sinh\lambda_n y + \delta_n y\sinh\lambda_n y\\ & {} + \int_{-1}^{1} k_n(y, \tilde y) c_n(\tilde y) d\tilde y\tag{22} \\ b_n = & A_n\cosh\mu_n y + B_n y\cosh\mu_n y + \Gamma_n\sinh\mu_n y + \Delta_n y\sinh\mu_n y\\ & {} + \int_{-1}^1 K_n(y, \tilde y) d_n(\tilde y) d\tilde y\tag{23} \end{align} for some functions $k$ and $K$. Coefficients $\alpha_n, \beta_n, \gamma_n$ and $\delta_n$ can be obtained by using the boundary conditions in the $y$-direction, and similarly for $A_n, B_n, \Gamma_n, \Delta_n$.
Now specialize to the case $g = 1$. We already have $c_n = \frac{1}{\left(n - \frac{1}{2}\right)\pi}$ and $d_n = 0$ from $(18)$ and $(19)$. (Because the non-homogeneous terms are constant, the method of undetermined coefficients is simpler than variation of parameters.) $c_n$ being constant makes $k_n = \frac{c_n}{2\lambda_n^4}$. Similarly, $d_n = 0$ makes $K_n = 0$. $a_n$ and $b_n$ become \begin{align} a_n = & \alpha_n\cosh\lambda_n y + \beta_n y\cosh\lambda_n y + \gamma_n\sinh\lambda_n y + \delta_n y\sinh\lambda_n y + \frac{c_n}{\lambda_n^4} \tag{S1}\\ b_n = & A_n\cosh\mu_n y + B_n y\cosh\mu_n y + \Gamma_n\sinh\mu_n y + \Delta_n y\sinh\mu_n y \tag{S2} \end{align} The boundary conditions in the $y$-direction have been translated to conditions for $a_n$ and $b_n$ in $(I1)$-$(I6)$.
All the conditions for $b_n$ are $0$, so $b_n = 0$.
For $a_n$, the conditions are $a_n(-1) = a_n(1) = 0$, $a''_n(-1) = c_n$ and $a''_n(1) = -c_n$. Plug the first two conditions into $(S1)$ an multiply by $\lambda_n^2$ to get \begin{align} -\frac{c_n}{\lambda_n^4} = & \alpha_n\cosh\lambda_n - \beta_n\cosh\lambda_n - \gamma_n\sinh\lambda_n + \delta_n\sinh\lambda_n \tag{C1}\\ -\frac{c_n}{\lambda_n^4} = & \alpha_n\cosh\lambda_n + \beta_n\cosh\lambda_n + \gamma_n\sinh\lambda_n + \delta_n\sinh\lambda_n \tag{C2} \end{align}
To use $a''_n(-1) = -c_n$ and $a''_n(1) = c_n$, compute $a''_n$ first: \begin{align} a'_n = & \alpha_n\lambda_n\sinh\lambda_n y + \beta_n \cosh\lambda_n y + \beta_n \lambda_n y\sinh\lambda_n y \\ & {} + \gamma_n\lambda_n\cosh\lambda_n y + \delta_n \sinh\lambda_n y + \delta_n\lambda_n y\cosh\lambda_n y \\ a''_n = & \alpha_n\lambda_n^2\cosh\lambda_n y + \beta_n\lambda_n\sinh\lambda_n y + \beta_n\lambda_n\sinh\lambda_n y + \beta_n\lambda_n^2 y\cosh\lambda_n y \\ & {} + \gamma_n\lambda_n^2\sinh\lambda_n y + \delta_n\lambda_n\cosh\lambda_n y + \delta_n\lambda_n\cosh\lambda_n y + \delta_n\lambda_n^2y\sinh\lambda_n y \\ = & \left(\alpha_n\lambda_n^2 + 2\delta_n\lambda_n\right)\cosh\lambda_n y + \beta_n\lambda_n^2 y\cosh\lambda_n y \\ & {} + \left(2\beta_n\lambda_n + \gamma_n\lambda_n^2\right)\sinh\lambda_n y + \delta_n\lambda_n^2y\sinh\lambda_n y \end{align} Then, plug in the values to get \begin{align} c_n = & \left(\alpha_n\lambda_n^2 + 2\delta_n\lambda_n\right)\cosh\lambda_n - \beta_n\lambda_n^2 \cosh\lambda_n \\ & {} - \left(2\beta_n\lambda_n + \gamma_n\lambda_n^2\right)\sinh\lambda_n + \delta_n\lambda_n^2\sinh\lambda_n \\ -c_n = & \left(\alpha_n\lambda_n^2 + 2\delta_n\lambda_n\right)\cosh\lambda_n + \beta_n\lambda_n^2 \cosh\lambda_n \\ & {} + \left(2\beta_n\lambda_n + \gamma_n\lambda_n^2\right)\sinh\lambda_n + \delta_n\lambda_n^2\sinh\lambda_n \end{align} Rearrange to get \begin{align} c_n = & \alpha_n \lambda_n^2 \cosh\lambda_n - \beta_n(\lambda_n^2 \cosh\lambda_n + 2\lambda_n\sinh\lambda_n)\\ & {} - \gamma_n \lambda_n^2 \sinh\lambda_n + \delta_n(\lambda_n^2\sinh\lambda_n + 2\lambda_n\cosh\lambda_n) \tag{C3}\\ -c_n = & \alpha_n\lambda_n^2\cosh\lambda_n + \beta_n(\lambda_n^2\cosh\lambda_n + 2\lambda_n\sinh\lambda_n)\\ & {} + \gamma_n\lambda_n^2\sinh\lambda_n + \delta_n(2\lambda_n\cosh\lambda_n + \lambda_n^2\sinh\lambda_n). \tag{C4} \end{align} Solve $(C1), (C2), (C3), (C4)$ for $\alpha_n, \beta_n, \gamma_n, \delta_n$. These equations can be simplified by dividing through by $\cosh\lambda_n$. This results in a system of linear equations \begin{align} \begin{pmatrix} 1 & -1 & -\tanh\lambda_n & \tanh\lambda_n \\ 1 & 1 & \tanh\lambda_n & \tanh\lambda_n \\ 1 & -1 - \frac{2}{\lambda_n}\tanh\lambda_n & -\tanh\lambda_n & \tanh\lambda_n + \frac{2}{\lambda_n} \\ 1 & 1 + \frac{2}{\lambda_n}\tanh\lambda_n & \tanh\lambda_n & \tanh\lambda_n + \frac{2}{\lambda_n}\end{pmatrix} \begin{pmatrix} \alpha_n\\ \beta_n\\ \gamma_n\\ \delta_n \end{pmatrix} = \frac{1}{\lambda_n^5\cosh\lambda_n} \begin{pmatrix} -1 \\ -1 \\ \lambda_n^2 \\ -\lambda_n^2 \end{pmatrix} \tag{SS} \end{align} Once $\alpha_n, \beta_n, \gamma_n$ and $\delta_n$ are known, put them back into $(S1)$ to get $a_n$. Go back to $(11)$ to get the solution $$ f(x, y) = \sum_{n=1}^\infty a_n(y) u_n(x). $$
Summary To compute the series solution, initialize the solution as zero. Then for each $n$,
Repeat for more accuracy.
Note: I am not 100% sure that all calculations are correct, but the method definitely is. If you spot any mistakes, please correct them.