# Solve inhomogeneous Poisson's equation, gives no result at all

I would like to solve the following two-dimensional inhomogeneous Poisson's equation in Mathematica including specific boundary conditions, and I know that an analytical solution exists, but Mathematica is not cooperating in this special case.

I would like to analytically solve the following Poisson's equation: $$\frac{\partial^2 u(x,y)}{\partial x^2}+\frac{\partial^2 u(x,y)}{\partial y^2}=6 x (1-y) y-2 x^3$$ with the boundary conditions $$u(0, y)=0,\ u(1, y)=(1-y)y,\ u(x, 0)=0,\ u(x, 1)=0$$ From a textbook I know that the solution is $$u(x,y)=y(1-y)x^3$$ If I try to solve this equation system by using DSolve

DSolve[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == -2 x^3 + 6 x (1 - y) y,
u[0, y] == 0,
u[1, y] == (1 - y) y,
u[x, 0] == 0,
u[x, 1] == 0}, u, {x, y}]


But Mathematica refused to do anything and is not returning any error message at all. Am I overlooking something or is such an PDE not solvable analytically in Mathematica?

Of course one could use NDSolve, but I would like to have the possibility to get more general solutions through DSolve.

-
This doc page says: DSolve can find the general solution for a restricted type of homogeneous linear second-order PDEs; namely, equations of the form a*D[u[x, y], {x, 2}] + b*D[u[x, y], x, y] + c*D[u[x, y], {y, 2}] == 0. Your equation has a nonvanishing non-principal part, and the algorithm used by DSolve is not applicable in this case. –  Sjoerd C. de Vries Jun 26 '13 at 20:13

DSolve is not equipped to solve this but it's pretty easy to solve. Given a polynomial $f$, you wish to find $u$ such that $$u_{xx}+u_{yy} = f(x,y),$$ together with polynomial boundary conditions on a square. It seems reasonably clear $u$ should be a multi-variate polynomial, so let's write down the general fifth degree multi-variate polynomial and interpret the PDE and boundary conditions as algebraic conditions on the coefficients.

Here are $f$ and $u$ defined in Mathematica.

Clear[f,u];
f[x_,y_] = 6x(1-y)y-2x^3;
u[x_,y_] = Sum[a[n,j]x^(n-j) y^j,{n,0,5},{j,0,n}]


Now, the PDE $u_{xx}+u_{yy} = f(x,y)$, can be expressed in terms of the coefficients. For example, the coefficient of $x^3$ in $f$ is $-2$. The same coefficient in $u_{xx}+u_{yy}$ can be found as follows:

Coefficient[D[u[x,y],x,x]+D[u[x,y],y,y],x^3]


Thus, we have $20a_{5,0}+2a_{5,2} = -2$. This is just one of a system of equations that must be satisfied. We can grab them all as follows.

rhs = Tuples[{Range[0,3],Range[0,3]}] /.
Append[CoefficientRules[6x(1-y)y-2x^3,{x,y}],{_,_} -> 0];
lhs = Tuples[{Range[0,3],Range[0,3]}] /.
Append[CoefficientRules[D[u[x,y],x,x]+D[u[x,y],y,y],{x,y}],{_,_} -> 0];


The Trues are a bit superfluous but don't hurt anything. Let's redefine $u$ by solving this system and plugging the resulting coefficients back into $u$.

u[x_,y_] = u[x,y] /. First[Solve[eqs]]


Note now that number of unknowns has been reduced - $a_{5,5}$, for example, appears three times. If this worked, then the following should be exactly $f(x,y)$.

D[u[x,y],x,x]+D[u[x,y],y,y]//Expand


Now, the solution can be narrowed down further using the boundary conditions. Consider, for example, $u(0,y)=0$:

u[0,y]


This is a polynomial in $y$ and can only equal zero if all the coefficients are zero. Thus, all the above coefficients must be zero. The same is true of $u(x,0)$ and $u(x,1)$. This quickly yields a huge reduction.

furtherReduction = #->0& /@ Cases[{u[0,y],u[x,0],u[x,1]},_a,Infinity]


If we plug these in, we should get even closer to the actual solution.

u[x,y] /. furtherReduction


In fact, we've already hit it exactly. Setting $x=1$, it's easy to see that the remaining boundary condition is satisfied.

-
perfect analysis. thank you for that. it shows again how well mathematica is suited in the academic study. (+1) –  Stefan Jun 26 '13 at 21:37
I'm impressed. Many thanks for this throughout analysis, I appreciate it! –  Rainer Jul 1 '13 at 18:15