# 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. – bbgodfrey 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. – Asatur Khurshudyan Jan 21 '17 at 20:57
• bbgodfrey, actually, the governing equation is nonlinear, so Fourier or other similar things are useful. – Asatur Khurshudyan Jan 21 '17 at 20:59
• You should add the b.c.s and P to your question by clicking the edit button. – xzczd 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: – Asatur Khurshudyan 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 – Asatur Khurshudyan 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} – Asatur Khurshudyan Jan 23 '17 at 7:13
• It doesnt work. – Asatur Khurshudyan Jan 23 '17 at 7:14
• @AsaturKhurshudyan Check my edit. – xzczd Jan 23 '17 at 8:21

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. – Asatur Khurshudyan Nov 26 '18 at 1:02
• What is here sol[x,y]`? – Alexei Boulbitch Feb 21 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. – user21 Feb 21 at 7:41