10
$\begingroup$

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?

$\endgroup$
8
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Jan 21, 2017 at 20:42
  • 1
    $\begingroup$ What are bc and P[x,y]? $\endgroup$
    – zhk
    Jan 21, 2017 at 20:52
  • 1
    $\begingroup$ 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. $\endgroup$ Jan 21, 2017 at 20:57
  • $\begingroup$ bbgodfrey, actually, the governing equation is nonlinear, so Fourier or other similar things are useful. $\endgroup$ Jan 21, 2017 at 20:59
  • 1
    $\begingroup$ You should add the b.c.s and P to your question by clicking the edit button. $\endgroup$
    – xzczd
    Jan 22, 2017 at 6:40

2 Answers 2

17
$\begingroup$

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,
  please find it in the link above.*)
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

Mathematica graphics

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,
  please find it in the link above. *)
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

Mathematica graphics

Notice setting proper initial values for FindRoot can be troublesome, but luckily it seems not to be a big problem in this case.

$\endgroup$
11
  • $\begingroup$ Actually, the deflection should be downwards. However, when I want to apply the code for the following set: $\endgroup$ Jan 23, 2017 at 7:13
  • $\begingroup$ 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 $\endgroup$ Jan 23, 2017 at 7:13
  • $\begingroup$ 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} $\endgroup$ Jan 23, 2017 at 7:13
  • 1
    $\begingroup$ Thanks, it is physically correct. Thank you very much. $\endgroup$ Jan 23, 2017 at 11:17
  • 1
    $\begingroup$ your answer solved my problem and, in fact, I am still using your code. So, there shouldn't be anything to be improved. $\endgroup$ Mar 20, 2020 at 10:07
15
$\begingroup$

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:

P[x_, y_] := x y
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}]
(* alternative workaround for version 13.0.1 *)
(* ufun = NDSolveValue[{eqn, bcs}, {u, v}, {x, 0, 1}, {y, 0, 1}][[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}]

enter image description here

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}]

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.