# Solving biharmonic equation with Mathematica

I would like to solve a biharmonic equation in polar coordinates of the form:

$\Delta \Delta \Phi[r,\theta] = r^i cos(j \, \theta) \quad i,j \, \epsilon \, \mathbb{N}_0$

I know that a solution for the homogeneous problem

$\Delta \Delta \Phi[r,\theta] = 0$

was found quite a while ago (100 years+), see here. I am wondering why Mathematica is not able to solve even the simplified case via

$DSolve[ \Delta[ \Delta \phi [r, \theta], \{r, \theta\}, "Polar"], \{r, \theta\}, "Polar"] == 0, \phi, \{r, \theta\}]$

Code to copy:

  DSolve[Simplify[
Laplacian[
Laplacian[ϕ[r, θ], {r, θ}, "Polar"], {r, θ}, "Polar"]] == 0, ϕ, {r, θ}]


because the problem is quite common in elasticity.

I know a similar question has been asked before, but the recommended solution is just valid for rotational symmetric case for which the PDE simplifies to an ODE.

I would be very thankful for any suggestion.

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Jun 22 '16 at 17:50
• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this meta Q&A helpful – Michael E2 Jun 22 '16 at 17:50

If you accept that the general solution can be constructed as an (infinite) Fourier series,

$$\sum_{m=-\infty}^{\infty}\phi_m(r) \exp(i m \theta),$$

then you can obtain the expression for $\phi_m(r)$ as follows:

eqn =
0 == Simplify@
Laplacian[ Laplacian[ ϕ[r] Exp[I m θ], {r, θ}, "Polar"], {r, θ}, "Polar"];

DSolve[eqn, ϕ, r]

(*
==> {{ϕ ->
Function[{r},
r^-m C[1] + r^(2 - m) C[2] + r^(2 + m) C[3] + r^m C[4]]}}
*)


The integration constants in this result can then be chosen differently for each $m$. This solves the homogeneous problem (if the constants are chosen to satisfy the given boundary conditions).

You can do the same for the inhomogeneous problem with specific $i, j$, e.g., $i=2$, $j=3$:

eqn1 =
r^2 Exp[I 3 θ] ==
Simplify@
Laplacian[Laplacian[ ϕ[r] Exp[I 3 θ], {r, θ}, "Polar"], {r, θ}, "Polar"];

DSolve[eqn1, ϕ, r]

(*
==> {{ϕ ->
Function[{r}, r^6/189 + C[1]/r^3 + C[2]/r + r^3 C[3] + r^5 C[4]]}}
*)

eqn2 =
r^2 Exp[-I 3 θ] ==
Simplify@
Laplacian[ Laplacian[ ϕ[r] Exp[-I 3 θ], {r, θ}, "Polar"], {r, θ}, "Polar"];

DSolve[eqn2, ϕ, r]

(*
==> {{ϕ ->
Function[{r}, r^6/189 + C[1]/r^3 + C[2]/r + r^3 C[3] + r^5 C[4]]}}
*)


Here, I broke the $\cos(3\theta)$ up into two parts containing $\exp(\pm 3 i\theta)$ and solved the two separately. Then you have to add both solutions and divide by two to get the solution to with $r^2 \cos(3 \theta)$ as the inhomogeneous term.

Finally, using the same superposition principle as above for some special choices of $i,j$, I tried the general solution for arbitrary $i,j$ in the following form:

eqn3 =
r^i Cos[j θ] ==
Simplify@
Laplacian[ Laplacian[ ϕ[r] Cos[j θ], {r, θ}, "Polar"], {r, θ}, "Polar"];

DSolve[eqn3, ϕ, r]

(*
==> {{ϕ ->
Function[{r},
r^(4 + i)/((-4 - i + j) (-2 - i + j) (2 + i + j) (4 + i + j)) +
r^-j C[1] + r^(2 - j) C[2] + r^(2 + j) C[3] + r^j C[4]]}}
*)


So for the inhomogeneous equation, all we had to do is choose the same $\theta$ dependence for the unknown solution as appears in the given right-hand side.

• Thank you so much, Jens! That's exactly what I was looking for and it works perfectly fine. – Paul Saturday Jun 24 '16 at 6:42