# 2D inhomogeneous biharmonic equation with wedged edge

I'm solving bending of rectangular plate while, boundary conditions are I have found similar problem solved: datavoreconsulting.com/programming-tips/numerically-solving-pdes-mathematica-finite-differences/ and used this method to solve my problem

Define Functions Make the domain into a grid

\[Nu] = 0.3;
ee = 210*10^9;
p = -19000
h = 7;
DD = (ee h^3)/(12 (1 - \[Nu]^2));
a = 1000;
b = 2 a;
xmin = -(1/2) a;
xmax = 1/2 a;
ymin = -(1/2) b;
ymax = 1/2 b;
xdivisions = 10;
ydivisions = 10;
dx = (xmax - xmin)/xdivisions
dy = (ymax - ymin)/ydivisions
q0 = p;
xgrid = Range[xmin, xmax, dx];
ygrid = Range[ymin, ymax, dy];
grid = Outer[{#1, #2} &, xgrid, ygrid];


Make an array of values for the solution. Each entry in the array corresponds to the value of the unknown function at a point in space

W = Array[w, {xdivisions + 1, ydivisions + 1}, {{xmin, xmax}, {ymin, ymax}}];

-19000

100

200


Specify boundary condition W=0

wl[y_] = 0;
wr[y_] = 0;
wb[x_] = 0;
wt[x_] = 0;


Find finite difference appoximations for different derivatives of u

dwdx = NDSolveFiniteDifferenceDerivative[{1, 0}, {xgrid, ygrid}, W];
dwdx2 = NDSolveFiniteDifferenceDerivative[{2, 0}, {xgrid, ygrid}, W];
dwdy = NDSolveFiniteDifferenceDerivative[{0, 1}, {xgrid, ygrid}, W];
dwdy2 = NDSolveFiniteDifferenceDerivative[{0, 2}, {xgrid, ygrid}, W];
dwdxdy = NDSolveFiniteDifferenceDerivative[{1, 1}, {xgrid, ygrid}, W];
dwdx4 = NDSolveFiniteDifferenceDerivative[{4, 0}, {xgrid, ygrid}, W];
dwdy4 = NDSolveFiniteDifferenceDerivative[{0, 4}, {xgrid, ygrid}, W];
dwdx2dy2 = NDSolveFiniteDifferenceDerivative[{2, 2}, {xgrid, ygrid}, W];


Boundary Conditions

(*Consider the left-side bc *)
leftbc = Table[dwdx[[1, k + 1]], {k, 0, Length[ygrid] - 1}];
leftbcw = Table[W[[1, k + 1]], {k, 0, Length[ygrid] - 1}];

(*Do the same for the right-side bc*)
rightbc = Table[dwdx[[-1, k + 1]], {k, 0, Length[ygrid] - 1}];
rightbcw = Table[W[[-1, k + 1]], {k, 0, Length[ygrid] - 1}];

(*And for the bottom side boundary condition*)
bottombc = Table[dwdy[[k + 1, 1]], {k, 0, Length[xgrid] - 1}];
bottombcw = Table[W[[k + 1, 1]], {k, 0, Length[xgrid] - 1}];

(*And for the top side bc*)
topbc = Table[dwdy2[[k + 1, -1]], {k, 0, Length[xgrid] - 1}];
topbcw = Table[W[[k + 1, -1]], {k, 0, Length[xgrid] - 1}];


Takes the left boundary conditions, and solves them to yield the values of \ u along the left side

   wleft = NSolve[Map[# == 0 &, leftbc], Table[W[[1, k]], {k, 1, Length[ygrid]}]];
wleftw = NSolve[Map[# == 0 &, leftbcw],
Table[W[[1, k]], {k, 1, Length[ygrid]}]];


takes the right boundary conditions, and solves them to yield the values of \ u along the right side

wright = NSolve[Map[# == 0 &, rightbc],
Table[W[[-1, k]], {k, 1, Length[ygrid]}]];
wrightw = NSolve[Map[# == 0 &, rightbcw],
Table[W[[-1, k]], {k, 1, Length[ygrid]}]];


takes the bottom boundary conditions, and solves them to yield the values \ of u along the bottom side

wbottom = NSolve[Map[# == 0 &, bottombc],
Table[W[[k, 1]], {k, 1, Length[xgrid]}]];
wbottomw =
NSolve[Map[# == 0 &, bottombcw], Table[W[[k, 1]], {k, 1, Length[xgrid]}]];


takes the top boundary conditions, and solves them to yield the values of u along the top side

 wtop = NSolve[Map[# == 0 &, topbc], Table[W[[k, -1]], {k, 1, Length[xgrid]}]];
wtopw = NSolve[Map[# == 0 &, topbcw],
Table[W[[k, -1]], {k, 1, Length[xgrid]}]];


To make sure boundary conditions are consitent, we only want one u value \ per point

wleft1 = wleft[[1, 1 ;; Length[ygrid]]];
wtop1 = wtop[[1, 2 ;; Length[xgrid] - 1]];
wright1 = wright[[1, 1 ;; Length[ygrid]]];
wbottom1 = wbottom[[1, 2 ;; Length[xgrid] - 1]];
wleft2 = wleftw[[1, 1 ;; Length[ygrid]]];
wtop2 = wtopw[[1, 2 ;; Length[xgrid] - 1]];
wright2 = wrightw[[1, 1 ;; Length[ygrid]]];
wbottom2 = wbottomw[[1, 2 ;; Length[xgrid] - 1]];
(*This is a list of all the boundary values of u*)
boundary1 = Join[wleft1, wright1, wtop1, wbottom1];
boundary2 = Join[wleft2, wright2, wtop2, wbottom2];

boundary = boundary2 /. boundary1;


Solving

We now create a set of equations using the PDE. We make a table, with each entry corresponding to an interior grid point. Each entry in the table becomes an equation, from the discretized PDE. We use our knowledge of the boundary conditions to eliminate the the values of u on the boundary

equations =
Map[(q0/DD == #) &,
Flatten[Table[
dwdx4[[i, j]] + 2*dwdx2dy2[[i, j]] + dwdy4[[i, j]], {i, 2,
Length[xgrid] - 1}, {j, 2, Length[ygrid] - 1}], 1] /. boundary /.
boundary2];

(*We solve these equations, for the interior values of u*)
intSol = NSolve[equations,
Flatten[W[[2 ;; Length[xgrid] - 1, 2 ;; Length[ygrid] - 1]]]][];


We substitute the interior values of u into the equations that determine the boundary values of u

boundarySol = boundary2 /. intSol;


solutionArray is an array of the values of u. I'm not sure why I had to do the last replacement rule twice!

solutionArray = ((W /. intSol) /. boundarySol) /. boundarySol;


Use this to make a table of 3d coordinates

dataPoints =
Table[{xmin + i*dx, ymin + j*dy, solutionArray[[i + 1, j + 1]]}, {i, 0,
Length[xgrid] - 1}, {j, 0, Length[ygrid] - 1}];
Min[solutionArray]

-90.0146


Plot the solution

ListPlot3D[Flatten[dataPoints, 1],
AxesLabel -> {Style[x, Medium, Blue], Style[y, Medium, Blue],
Style[u, Medium, Blue]}, PlotRange -> All]


but solution differs from FEM softvare solution and some other mathematica solution, my code gives max deflection 90.0146 but it should be around 1.5. Am I applying boundary conditions wrong way or does mathematica method give huge error?

I have solved the problem by programming my own finite difference parameters.

new finite difference approximation functions:
dwdx = Table[(W[[i + 1, j]] - W[[i - 1, j]])/(2*dx), {i, 3,
Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdy = Table[(W[[i, j + 1]] - W[[i, j - 1]])/(2*dy), {i, 3,
Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdx2 = Table[(W[[i + 1, j]] - 2 W[[i, j]] +
W[[i - 1, j]])/(dx^2), {i, 3, Length[xgrid] - 2}, {j, 3,
Length[ygrid] - 2}];
dwdy2 = Table[(W[[i, j + 1]] - 2 W[[i, j]] +
W[[i, j - 1]])/(dy^2), {i, 3, Length[xgrid] - 2}, {j, 3,
Length[ygrid] - 2}];
dwdx4 = Table[(W[[i + 2, j]] - 4 W[[i + 1, j]] + 6 W[[i, j]] -
4 W[[i - 1, j]] + W[[i - 2, j]])/(dx^4), {i, 3,
Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdy4 = Table[(W[[i, j + 2]] - 4 W[[i, j + 1]] + 6 W[[i, j]] -
4 W[[i, j - 1]] + W[[i, j - 2]])/(dy^4), {i, 3,
Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdx2dy2 =
Table[(W[[i + 1, j + 1]] - 2 W[[i + 1, j]] + W[[i + 1, j - 1]] -
2 W[[i, j + 1]] + 4 W[[i, j]] - 2 W[[i, j - 1]] +
W[[i - 1, j + 1]] - 2 W[[i - 1, j]] + W[[i - 1, j - 1]])/(dx^2*
dy^2), {i, 3, Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];


Now I have another problem to solve, bending of a plate with hole

• I've solved a very similar problem based on pdetoae (Which is a shell of FiniteDifferenceDerivative) here, you can have a look. – xzczd Sep 15 '17 at 14:49
• Thnx i'll take a look when i get to my comp. – Katarina Sep 15 '17 at 17:49
• @xzczd I have solved my problem but new problem has occured, and I do not understand your solution. can you solve the same problem as in my question, using your algorithm so I can compare results? – Katarina Oct 20 '17 at 16:16

## 1 Answer

I haven't check your code, but here's my solution based on my pdetoae. You may feel it confusing at first glance, but I believe you'll be able to understand it if you settle down and read through it. The idea is simple: generate the difference equations automatically with ptoafunc, then remove those redundant equations with del, extract the coefficients of the linear equations system with CoefficientArrays and solve it with LinearSolve (you can also solve the equation system with Solve directly, but it'll be much slower), and built a InterpolatingFunction with ListInterpolation.

ν = 3/10;
Ε = 210*10^9;
qz = -19000;
h = 7;
\[ScriptCapitalD] = (Ε h^3)/(12 (1 - ν^2));
a = 1000;
b = 2 a;

lap = Laplacian[#, {x, y}] &;

With[{w = w[x, y]},
eq = lap@lap@w == qz/\[ScriptCapitalD];
(*bc={
{w == 0,D[w,x,x] == 0}/.{{x -> 0},{x -> a}},
{{w == 0,D[w,y,y] == 0}/.y -> 0,{w == 0,D[w,y] == 0}/.y -> b}}*)
bc = {
{w == 0, D[w, x] == 0} /. {{x -> 0}, {x -> a}},
{{w == 0, D[w, y] == 0} /. y -> 0, {w == 0, D[w, y, y] == 0} /. y -> b}};]

points = 50;
{grid@x, grid@y} = Array[# &, points, #] & /@ {{0, a}, {0, b}};

difforder = 4;
ptoafunc = pdetoae[w[x, y], grid /@ {x, y}, difforder];
del = #[[3 ;; -3]] &;

ae = del /@ del@ptoafunc@eq;
aebc = MapAt[del, ptoafunc@bc, {1, All, All}];

var = Outer[w, grid@x, grid@y] // Flatten;
{blst, mat} = CoefficientArrays[{ae, aebc} // Flatten, var];
sollst = LinearSolve[mat, -N@blst];
sol = ListInterpolation[Partition[sollst, points], {grid@x, grid@y}];

Plot3D[sol[x, y], {x, 0, a}, {y, 0, b}] Please notice the b.c.s in your first and last picture are different, and it's not clear to me which boundary is wedged in the last case so I chose one arbitrarily. Anyway the result agrees well with the one in the picture:

NMinimize[{sol[x, y], 0 < x < a, 0 < y < b}, {x, y}]
(* {-7.43643, {x -> 500., y -> 1086.09}} *)


# Update

If you feel confused about del, the following is an alternative method without it:

fullsys = ptoafunc@{eq, bc};
{blst, mat} = CoefficientArrays[fullsys // Flatten, var];
sollst = LeastSquares[mat, -N@blst, Method -> Direct]; // AbsoluteTiming
sol = ListInterpolation[Partition[sollst, points], {grid@x, grid@y}];

• xzczd@ thnx, I was solving two cases, first case: 3 freely supported, 1 wedged edge and second case: 1 freely supported, 3 wedged edge, I'll change picture as soon as possible. I had to give classes so my research was left aside, but hope to do some work in between my lectures. – Katarina Nov 15 '17 at 15:33
• I have changed boundary conditions and picture to correct picture – Katarina Nov 15 '17 at 15:53