# 2D inhomogeneous biharmonic equation

How to solve a 2D inhomogeneous biharmonic equation with NDSolve?

I tried the following code:

P[x_, y_] := x y
eq = Laplacian[Laplacian[w[x, y], {x, y}], {x, y}] == x*y;
bc = {w[0, y] == w[1, y] == w[x, 0] == w[x, 1] == 0,
Derivative[2, 0][w][0, y] == Derivative[2, 0][w][1, y] ==
Derivative[0, 2][w][x, 0] == Derivative[0, 2][w][x, 1] == 0};
NDSolve[{eq == P[x, y], bc}, w, {x, 0, 1}, {y, 0, 1}]


but it says

NDSolve::femcmsd: The spatial derivative order of the PDE may not exceed two.

How to derive the solution?

• As the error message says, NDSolve is not able to solve this problem as written. However, you could Fourier transform the system in one or both dimensions and proceed from there. Jan 21 '17 at 20:42
• What are bc and P[x,y]?
– zhk
Jan 21 '17 at 20:52
• MMM, bcs are w[0,y]=w[1,y]=w[x,0]=w[x,1]=0 and (D[w[x,y],x,x]/.x->0)=(D[w[x,y],x,x]/.x->1)=(D[w[x,y],y,y]/.y->0)=(D[w[x,y],y,y]/.y->1)=0 and $P[x,y]=x*y$, for example. Jan 21 '17 at 20:57
• bbgodfrey, actually, the governing equation is nonlinear, so Fourier or other similar things are useful. Jan 21 '17 at 20:59
• You should add the b.c.s and P to your question by clicking the edit button. Jan 22 '17 at 6:40

As mentioned in the warning, currently "FiniteElement" method can't handle 4th order spatial derivatives. So let me show you a FDM-based solution. I'll use pdetoae for the generation of difference equation:

P[x_, y_] := x y
eq = Laplacian[Laplacian[w[x, y], {x, y}], {x, y}] == P[x, y];
bc = {w[0, y] == w[1, y] == w[x, 0] == w[x, 1] == 0,
Derivative[2, 0][w][0, y] == Derivative[2, 0][w][1, y] ==
Derivative[0, 2][w][x, 0] == Derivative[0, 2][w][x, 1] == 0} /.
Equal[a__, b_] :> Thread[{a} == b];
{bcy, bcx} = GatherBy[Flatten@bc, FreeQ[#, _[0 | 1, y]] &];
domain = {0, 1};
points = 25;
grid = Array[# &, points, domain];
difforder = 4;
(*Definition of pdetoae isn't included in this code piece,
ptoafunc = pdetoae[w[x, y], {grid, grid}, difforder];
var = Outer[w, grid, grid] // Flatten;

del = #[[3 ;; -3]] &;

ae = del /@ del@ptoafunc@eq;
aebcx = ptoafunc@bcx;
aebcy = del /@ ptoafunc@bcy;

{b, m} = CoefficientArrays[{ae, aebcx, aebcy} // Flatten, var];

sollst = LinearSolve[m, -N@b];


Remark

If you have difficulty in understanding the usage of del, the following is an alternative way for calculating sollst:

fullsys = ptoafunc@{eq, bcx, bcy} // Flatten;
{b, m} = CoefficientArrays[fullsys, var];
sollst = LeastSquares[m, -N@b]; // AbsoluteTiming


Notice this approach is slower.

sol = ListInterpolation[Partition[sollst, points], {grid, grid}];

Plot3D[sol[x, y], {x, ##}, {y, ##}] & @@ domain Notice I've modified the definition of bc because pdetoae can't parse continued equality at the moment i.e. something like a == b == c isn't supported yet.

Solution for the problem in the comments below

The new-added example in the comment has a nonlinear inhomogeneous term, so LinearSolve can't be used any more, we can use FindRoot instead:

nu = 0.33; h = 0.01; Ye = 2 10^11; P1 = 10^5;
N11[x_, y_] = (Ye h)/(2 (1 - nu^2)) ((D[w[x, y], x])^2 + nu (D[w[x, y], y])^2);
N22[x_, y_] = (Ye h)/(2 (1 - nu^2)) (nu (D[w[x, y], x])^2 + (D[w[x, y], y])^2);
N12[x_, y_] = (Ye h)/(2 (1 + nu)) D[w[x, y], x] D[w[x, y], y] ;
P[x_, y_] =
N11[x, y] D[w[x, y], x, x] - N22[x, y] D[w[x, y], y, y] -
2 N12[x, y] D[w[x, y], x, y] - P1;
eq = (Ye h^3)/(12 (1 - nu^2)) Laplacian[Laplacian[w[x, y], {x, y}], {x, y}] == -P[x,
y]; bc = {w[x, 0] == w[x, 1] == 0,
Derivative[2, 0][w][0, y] == Derivative[2, 0][w][1, y] == 0,
Derivative[0, 2][w][x, 0] == Derivative[0, 2][w][x, 1] ==
0, (Ye h^3)/(12 (1 - nu^2)) (Derivative[3, 0][w][0, y] +
2 Derivative[1, 2][w][0, y]) + P1 Derivative[1, 0][w][0, y] ==
0, (Ye h^3)/(12 (1 - nu^2)) (Derivative[3, 0][w][1, y] +
2 Derivative[1, 2][w][1, y]) + P1 Derivative[1, 0][w][1, y] == 0} /.
Equal[a__, b_] :> Thread[{a} == b];
{bcy, bcx} = GatherBy[Flatten@bc, FreeQ[#, _[0 | 1, y]] &];
domain = {0, 1};
points = 25;
grid = Array[# &, points, domain];
difforder = 4;
(* Definition of pdetoae isn't included in this code piece,
ptoafunc = pdetoae[w[x, y], {grid, grid}, difforder];
del = #[[3 ;; -3]] &;
ae = del /@ del@ptoafunc@eq;
aebcx = ptoafunc@bcx;
aebcy = del /@ ptoafunc@bcy;
var = Outer[w, grid, grid] // Flatten;

solrule = FindRoot[Rationalize[{ae, aebcx, aebcy} // Flatten, 0], {#, 0} & /@ var,
WorkingPrecision -> 16]; // AbsoluteTiming
sollst = Replace[solrule, (w[x_, y_] -> z_) :> {x, y, z}, {1}];
sol = Interpolation@sollst;
Plot3D[sol[x, y], {x, ##}, {y, ##}] & @@ domain Notice setting proper initial values for FindRoot can be troublesome, but luckily it seems not to be a big problem in this case.

• Actually, the deflection should be downwards. However, when I want to apply the code for the following set: Jan 23 '17 at 7:13
• nu = 0.33; h = 0.01; Ye = 2 10^11; P1 = 10^5; N11[x_, y_] := (Ye h)/( 2 (1 - nu^2)) ((D[w[x, y], x])^2 + nu (D[w[x, y], y])^2) N22[x_, y_] := (Ye h)/( 2 (1 - nu^2)) (nu (D[w[x, y], x])^2 + (D[w[x, y], y])^2) N12[x_, y_] := (Ye h)/(2 (1 + nu)) D[w[x, y], x] D[w[x, y], y] P[x_, y_] := N11[x, y] D[w[x, y], x, x] - N22[x, y] D[w[x, y], y, y] - 2 N12[x, y] D[w[x, y], x, y] - P1 Jan 23 '17 at 7:13
• eq = (Ye h^3)/(12 (1 - nu^2)) Laplacian[Laplacian[w[x, y], {x, y}], {x, y}] == -P[x, y]; bc = {w[x, 0] == w[x, 1] == 0, Derivative[2, 0][w][0, y] == Derivative[2, 0][w][1, y] == 0 == Derivative[0, 2][w][x, 0] == Derivative[0, 2][w][x, 1] == 0, (Ye h^3)/( 12 (1 - nu^2)) (Derivative[3, 0][w][0, y] + 2 Derivative[1, 2][w][0, y]) + P1 Derivative[1, 0][w][0, y] == 0, (Ye h^3)/( 12 (1 - nu^2)) (Derivative[3, 0][w][1, y] + 2 Derivative[1, 2][w][1, y]) + P1 Derivative[1, 0][w][1, y] == 0} Jan 23 '17 at 7:13
• Thanks, it is physically correct. Thank you very much. Jan 23 '17 at 11:17
• your answer solved my problem and, in fact, I am still using your code. So, there shouldn't be anything to be improved. Mar 20 '20 at 10:07

Update:

The example has been added to the help system. You can find it by clicking on the message NDSolve::femcmsd and following the link or by going to FEMDocumentation/ref/message/InitializePDECoefficients/femcmsd

For completeness I'd like to show that you can use the FEM to solve the biharmonic equation. The trick is to rewrite the 4th order equation as a system of two second order equations like so:

eqn = {Laplacian[u[x, y], {x, y}] == v[x, y],
Laplacian[v[x, y], {x, y}] == P[x, y]};
bcs = {u[0, y] == u[1, y] == u[x, 0] == u[x, 1] == 0,
v[0, y] == v[1, y] == v[x, 0] == v[x, 1] == 0};
ufun = NDSolveValue[{eqn, bcs}, u, {x, 0, 1}, {y, 0, 1}]


Note that the derivative boundary conditions from the original problem now are dirichlet conditions for the system of the equations.

A plot and a comparison to the other solutions show that is works well:

Plot3D[ufun[x, y], {x, 0, 1}, {y, 0, 1}] Compare this answer (ufun) to the answer given in xzczd's post (sol) to show that they match up.

Plot3D[ufun[x, y] - sol[x, y], {x, 0, 1}, {y, 0, 1}] • Really nice and simpler. Nov 26 '18 at 1:02
• What is here sol[x,y]? Feb 21 '19 at 7:31
• @AlexeiBoulbitch, the solution from the other answer. To sow that the answers match up. I'll edit the post to make that clearer. Feb 21 '19 at 7:41
• @use21 can we solve 4th order problems using Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> "TensorProductGrid"}}? Jul 20 '20 at 16:40
• @ABCDEMMM, time dependent yes, stationary like this one, no. Jul 21 '20 at 4:58