8
$\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$
  • $\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 '17 at 20:42
  • $\begingroup$ What are bc and P[x,y]? $\endgroup$ – zhk Jan 21 '17 at 20:52
  • $\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$ – Asatur Khurshudyan Jan 21 '17 at 20:57
  • $\begingroup$ bbgodfrey, actually, the governing equation is nonlinear, so Fourier or other similar things are useful. $\endgroup$ – Asatur Khurshudyan Jan 21 '17 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 '17 at 6:40
15
$\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$
  • $\begingroup$ Actually, the deflection should be downwards. However, when I want to apply the code for the following set: $\endgroup$ – Asatur Khurshudyan Jan 23 '17 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$ – Asatur Khurshudyan Jan 23 '17 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$ – Asatur Khurshudyan Jan 23 '17 at 7:13
  • $\begingroup$ It doesn`t work. $\endgroup$ – Asatur Khurshudyan Jan 23 '17 at 7:14
  • $\begingroup$ @AsaturKhurshudyan Check my edit. $\endgroup$ – xzczd Jan 23 '17 at 8:21
10
$\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:

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

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$
  • 1
    $\begingroup$ Really nice and simpler. $\endgroup$ – Asatur Khurshudyan Nov 26 '18 at 1:02
  • $\begingroup$ What is here sol[x,y]? $\endgroup$ – Alexei Boulbitch Feb 21 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 at 7:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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