2D inhomogeneous biharmonic equation

Question

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
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 06:40
Possible duplicate of Dynamic Euler–Bernoulli beam equation — user21, Jan 23 '17 at 08:46
I'm voting to reopen this question because the 4 votes asking for additional information is actually voted before OP adds specific code i.e. additional information. (Notice the decisive vote by @user21 is suggesting a possible duplicate. ) — xzczd, Aug 11 '18 at 15:15

xzczd · Accepted Answer · 2019-02-21T03:34:05.063

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.

Actually, the deflection should be downwards. However, when I want to apply the code for the following set: — Asatur Khurshudyan, Jan 23 '17 at 07: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 07: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 07:13
To the downvoter, I am interested in what's wrong with my answer, would you please elaborate. I'm not trying to complain here, I just want to improve my answer if possible. — xzczd, Mar 19 '20 at 09:03
your answer solved my problem and, in fact, I am still using your code. So, there shouldn't be anything to be improved. — Asatur Khurshudyan, Mar 20 '20 at 10:07
@xzczd can we use FDM-based solution for 2D/3D solid mechanic problems in MMA? If Yes, I will give a try in MMA! — ABCDEMMM, Jun 20 '20 at 14:33
@ABCDEMMM As long as the region is regular. (Rectangle, ball, etc. ) Somewhat related: https://scicomp.stackexchange.com/q/10187/5331 — xzczd, Jun 20 '20 at 15:02

user21 · Answer 2 · 2022-05-04T04:46:38.067

15

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

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

edited May 04 '22 at 04:46

answered Nov 07 '18 at 14:09

user21

39,710
8
110
167

3

Really nice and simpler. – Asatur Khurshudyan Nov 26 '18 at 01:02
1

What is here sol[x,y]? – Alexei Boulbitch Feb 21 '19 at 07: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 '19 at 07:41
@use21 can we solve 4th order problems using Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> "TensorProductGrid"}}? – ABCDEMMM Jul 20 '20 at 16:40
@ABCDEMMM, time dependent yes, stationary like this one, no. – user21 Jul 21 '20 at 04:58
@user21 thanks! If I want to solve a 4th order One-Dimensional problems, is it possible direct using NDSolve ? – ABCDEMMM Jul 21 '20 at 08:14
@ABCDEMMM, depending on the problem, depending on the definition of 'direct', yes. – user21 Jul 21 '20 at 08:18
@user21 is it that NSolve can only solve second order ode problems? example, plate or beam? – ABCDEMMM Jul 21 '20 at 08:24
@ABCDEMMM, I assume you mean NDSolve not NSolve. Yes, NDSolve can solve other then second order problems. The documentation has many examples. I do not understand what you want to know with: 'example, plate or beam ?' – user21 Jul 21 '20 at 08:46
@ABCDEMMM, do you mean this – user21 Jul 21 '20 at 08:49
@user21 yes, I am thinking, can we use mixed formulation for this problem, two two second gradient fields. – ABCDEMMM Jul 21 '20 at 08:51
@ABCDEMMM, unfortunately I do not understand what you are asking. – user21 Jul 21 '20 at 08:52
@user21 I cannot seem to run the first block of code in the answer to solve for u[x,y] in Mathematica 13.0.1, as in it does not evaluate to solve for u. Has something changed ? – Charlie May 03 '22 at 20:00
@Charlie, try to solve for u and v and not only for u. – user21 May 03 '22 at 20:12
@user21 I gave NDSolve the arguments {u,v} in sol = NDSolveValue[{eqn, bcs}, {u,v}, {x, 0, 1}, {y, 0, 1}] and it does give a solution but then ufun does not match sol[[1]] – Charlie May 03 '22 at 20:18
@Charlie, works just fine for me. – user21 May 04 '22 at 04:44

2D inhomogeneous biharmonic equation

2 Answers2

Linked

Related