# Coupled PDE using NDSolve

The problem is to study the solution of a biharmonic which can be decomposed into a coupled PDE system $$\Delta f - \nabla f\cdot\Gamma_{1} = g \\ \Delta g = \nabla f\cdot\Gamma_{2}$$where $$f, g$$ are smooth in $$\Omega = [-1, 1]\times[-1, 1]$$, and $$\Gamma_ 1$$ and $$\Gamma_ 2$$ are continuous in $$\Omega$$ and
$$\Gamma_ 1 = \frac{1}{(x^2+y^2)^{3/2}}(x, y)$$ and $$\Gamma_ 2 = \frac{1}{(x^2+y^2)^{3/2}}(-y, x)$$ .The boundary conditions are : $$\psi (-1, y) = 0, \quad\psi (1, y) = 1, \quad\hat {n}\cdot \nabla\psi | _ {y = \pm 1} = 0$$ I am not sure if I get the right solution because it gives me $$f = x - y$$ and $$g = x + y$$ which does not satisfy the PDE system. I don' t trust the solutions I get, because they remain unchanged when I change the vectors $$\Gamma$$ .

eqn1 = {Laplacian[psi[x, y], {x, y}] -
1/(x^2 +
y^2) (D[psi[x, y], {x, 1}] x/Sqrt[x^2 + y^2] +
D[psi[x, y], {y, 1}] y/Sqrt[x^2 + y^2]) == u1[x, y],
Laplacian[u1[x, y], {x, y}] - 1/(
x^2 + y^2) (D[psi[x, y], {x, 1}] y/Sqrt[x^2 + y^2] -
D[psi[x, y], {y, 1}] x/Sqrt[x^2 + y^2]) == 0 };
bcs = {psi[-1, y] == 1, psi[1, y] == 0}
sol = NDSolveValue[{eqn1, bcs}, {u1, psi}, {x, -1, 1}, {y, -1, 1}]


Exeprts, I am not in need of the full solution or anything. I am just in need of a hint as to whether this is the right approach to solve this problem, or is finite difference the way to go about solving a fourth order PDE. Thanks for the help!

EDIT :: (Apologies for editing this question, instead of starting a new one if that's the policy.)

I am confused on how to implement boundary conditions on a coupled system of PDEs. Let me state the physical problem that I am solving which will help explain my boundary conditions for the problem. The PDE system stated in the first equation describes the hydrodynamic flow where $$f$$ is the stream function and the fluid velocity is related to the stream function $$f$$ by $$v = \hat{z}\times \nabla f$$.This is done to impose the incompressibility condition that $$\nabla\cdot v=0$$ My boundary conditions are that the fluid enters from the right and exits from the left, and that on the top and bottom of the box the fluid does not slip, so that is the tangential component of $$v$$ is 0.

My idea was that I can set a Dirichlet condition on the top and bottom of the box so that the equipotential lines are horizontal in the direction of the fluid flow, and on the left and right hand side I can set a uniform flux in and out parallel to the horizontal axis, as shown in the figure in 1. How do I go about setting a Neumann condition on $$u_1$$. If I set the second equation equal to a Neumann value Mathematica ignores it and tells me that no boundary conditions have been set. Perhaps then I would need to set a Dirichlet condition on the left and right also, that is $$f(x=\pm L_x/2,y) = c_0 y + c_1$$? A hint in the right direction would be enough. Thanks for the help.

EDIT2 : Poiseuille flow :

T = 500; \[Nu] = 0.1 ;  t0 = 0.01;

\[Psi]1[0][x_, y_] := 0 ;

Do[{\[Psi]1[i], \[Xi]1[i]} =
NDSolveValue[{{-(\[Psi][x, y] - \[Psi]1[i - 1][x, y])/
t0 + \[Nu] Laplacian[\[Psi][x, y], {x, y}] - \[Xi][x, y],
Laplacian[\[Xi][x, y], {x, y}]} == {0,
0}, {DirichletCondition[\[Psi][x, y] == 1,
y == Ly/2 || -Lx/2 <= x <= Lx/2],
DirichletCondition[\[Psi][x, y] == -1,
y == -Ly/2 || -Lx/2 <= x <= Lx/2],
DirichletCondition[\[Xi][x, y] == 0, x == Lx/2],
DirichletCondition[\[Xi][x, y] == 0, x == -Lx/2],
DirichletCondition[\[Psi][x, y] == 2 y/Ly,
x == Lx/2]}}, {\[Psi], \[Xi]}, {x,
y} \[Element] \[CapitalOmega],
Method -> {"FiniteElement",
"InterpolationOrder" -> {\[Psi] -> 2, \[Xi] -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}], {i, 1, T}];


How to take care of $$\xi$$ ?

• In Latex you define solution on $\Omega$=Rectangle[{-1,-1},{1,1}], while in NDSolve you use one part $\Omega$ only {x, 0, 1}, {y, 0, 1}. Is it typo? May 4 at 2:14
• @AlexTrounev yes sorry for the typo. May 4 at 6:13
• I don't understand how you model has been derived from NSE. Could you give a link to the paper where this model described? May 7 at 5:24
• here's a reference - escholarship.org/uc/item/8kw631r4 May 7 at 8:03
• Ah, this is MHD, and therefore part of this system describes magnetic field. Could you show what part of your model follows from equations (2), (3) of the paper linked? May 7 at 9:36

To solve this problem we can use FEM as follows

Clear["Global*"]
Needs["NDSolveFEM"]
reg = Rectangle[{-1, -1}, {1, 1}]; mesh =
ToElementMesh[reg, MaxCellMeasure -> 0.001]

eqn1 = {Laplacian[psi[x, y], {x, y}] -
1/(x^2 + y^2) (D[psi[x, y], {x, 1}] x/Sqrt[x^2 + y^2] +
D[psi[x, y], {y, 1}] y/Sqrt[x^2 + y^2]) == u1[x, y],
Laplacian[u1[x, y], {x, y}] -
1/(x^2 + y^2) (D[psi[x, y], {x, 1}] y/Sqrt[x^2 + y^2] -
D[psi[x, y], {y, 1}] x/Sqrt[x^2 + y^2]) == 0};
bcs = {DirichletCondition[psi[x, y] == 0, x == -1],
DirichletCondition[psi[x, y] == 1, x == 1]};
sol = NDSolve[{eqn1, bcs}, {u1, psi}, Element[{x, y}, mesh]]


Since there is no any boundary condition for u1 we have a message

NDSolve::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {u1}; the result may not be unique.


Visualization

{Plot3D[Evaluate[psi[x, y] /. sol[[1]]], Element[{x, y}, mesh],
ColorFunction -> "Rainbow", MeshStyle -> White,
PlotTheme -> "Marketing", Boxed -> False, AxesLabel -> Automatic,
PlotLabel -> "psi"],
Plot3D[Evaluate[u1[x, y] /. sol[[1]]], Element[{x, y}, mesh],
ColorFunction -> "Rainbow", MeshStyle -> White,
PlotTheme -> "Marketing", Boxed -> False, AxesLabel -> Automatic,
PlotLabel -> "u1"]}


Update 1. In a case of homogenous flow on the right and left borders we have

Clear["Global*"]
Needs["NDSolveFEM"]
reg = Rectangle[{-1, -1}, {1, 1}]; mesh =
ToElementMesh[reg, MaxCellMeasure -> 0.001];

eqn1 = {Laplacian[psi[x, y], {x, y}] -
1/(x^2 + y^2) (D[psi[x, y], {x, 1}] x/Sqrt[x^2 + y^2] +
D[psi[x, y], {y, 1}] y/Sqrt[x^2 + y^2]) - u1[x, y] == 0,
Laplacian[u1[x, y], {x, y}] -
1/(x^2 + y^2) (D[psi[x, y], {x, 1}] y/Sqrt[x^2 + y^2] -
D[psi[x, y], {y, 1}] x/Sqrt[x^2 + y^2]) == 0};
bcs = {DirichletCondition[{psi[x, y] == y,
u1[x, y] == -y/(x^2 + y^2)^(3/2)}, y == -1 || y == 1],
DirichletCondition[psi[x, y] == y, x == 1 || x == -1]};
sol = NDSolve[{eqn1, bcs}, {u1, psi}, Element[{x, y}, mesh]]


Visualization

StreamDensityPlot[
Evaluate[{-Derivative[0, 1][psi][x, y],
Derivative[1, 0][psi][x, y]} /. sol[[1]]], Element[{x, y}, mesh],
ColorFunction -> Hue, PlotLegends -> Automatic]


• Thanks for the answer Alex. I had a question related to this problem, when imposing boundary conditions for a coupled system, is it possible to define a Dirichlet and Neumann boundary condition simultaneously for the two dependent variables in the system. What I mean is that in the problem above the psi is the stream function that imposes an incoming and outgoing stream flux on the top and bottom of the box and u1 is the laplacian of this quantity which is zero, so is it possible to impose a Dirichlet boundary condition of u1=0 and a Neumann boundary condition of constant on psi ? May 6 at 19:53
• @Charlie Do you mean condition for velocity $\vec{v}=0$ at $y=\pm 1$? May 9 at 10:53
• @ Alex, the boundary conditions on $y=\pm 1$ are the no-slip conditions May 9 at 15:54
• @Charlie It is clear from Figure 1 you shown that you mean slip conditions $\vec{v}.\vec{n}=0$. For the no-slip conditions in standard approach we need to add $f=0, \partial^2 f/\partial n^2=0$ on the border. But your equations are not in standard form. Can you show your no-slip boundary conditions in mathematical form? Above you suppose $f=0, u1=0$, but it is not working as no-slip conditions. May 10 at 3:41
• apologies I didn't specify the boundary conditions properly, I think the boundary conditions should be - at the top and bottom $\psi=1$ and $\psi=0$ respectively, which implies $v=\hat{z}\times\nabla\psi = 0$ and thus $u_1 = 0$which imposes the no-slip condition. At the left and right, flow is homogenous and $\psi = y$ and $u_1$ is found from second PDE in the coupled PDE system. May 10 at 16:46