I was trying to solve a nonlinear PDE

$-\Delta u(x,y) = u(x,y)^2$

in the domain $[0,1] \times[0,1]$ with Dirichlet boundary condition $u(x,y) =1$ on the boundary.

The linearized PDE is

$\nabla\cdot(-\nabla du-\nabla u^k) = (u^k+du)^2 \approx (u^k)^2+2u^k du.$

Rearrange the above equation into standard form according to the details of InitializePDECoefficients in the Mathematica 11.0 documentation,

$\nabla\cdot(-\nabla du-\nabla u^k) + (-2u^k) du - (u^k)^2 =0.$

According to the details of InitializePDECoefficients, enter image description here

the corresponding Coefficients are:

LoadCoefficients $f_1 \to (u^k)^2$

ReactionCoefficients $a_{11} \to -2u^k$

LoadDerivativeCoefficients $\gamma_1 \to \nabla u^k \, (i.e.,\{\frac{\partial u^k}{\partial x},\frac{\partial u^k}{\partial y}\})$.

The numerical solution blows up if I use the above coefficients. It will converge if I change the sign of the LoadDerivativeCoefficients to:

LoadDerivativeCoefficients $\gamma_1 \to -\nabla u^k \, (i.e.,\{-\frac{\partial u^k}{\partial x},-\frac{\partial u^k}{\partial y}\})$.

I'm confused why it works by changing the sign of the LoadDerivativeCoefficients. I also do not understand how the boundary condition is applied in the code. Can anyone explain it to me?


According to @user21, there is a typo in the ref page of InitializePDECoefficients. The sign of $\alpha$ and $\gamma$ should be switched.

My understanding of the second question according to @user21's ANSWER: To implement the linearized equation using FEM in Mathematica, we break it into two parts, that is, the linear part and the non-linear ($u^k$-dependent) part:

$\nabla\cdot(-\nabla du)$ = 0 -- (1)


$\nabla\cdot(-\nabla u^k) + (-2u^k) du - (u^k)^2 =0.$ -- (2)

We create the load vector and stiffness matrix of the (1) outside the Newton loop to avoid repeated calculations of this part. Inside the loop, we update the load vector and stiffness matrix of the (2) for previous solution $u^k$. The total load vector and stiffness matrix are obtained by adding up the two corresponding parts. Since the increment $du$ should satisfy the Dirichlet BC that $du=0$ on the boundary, we then need to modify the load vector and stiffness matrix accordingly. It is important that the initial guess $u^0$ should satisfy the Dirichlet BC that $u=0$ on the boundary.

  • $\begingroup$ Well, how have we learnt to solve nonlinear equations? Hm. Newton's method? $\endgroup$ – Henrik Schumacher Oct 23 '17 at 21:49
  • $\begingroup$ Thank you for the reply. @HenrikSchumacher. I understand the nonlinear loop part using Newton method in general, as I have denoted u^k as the previous solution in linearizing the nonlinear PDE. However, I'm not clear about the two points: (1) Why the sign of LoadDerivativeCoefficients is negative? and (2) Why the Dirichlet boundary condition is applied in the linearized part? Should the boundary condition of $du$ be 0? $\endgroup$ – Wilhelm Oct 23 '17 at 22:04

That is an unfortunate typo in the documentation in V11.2. The correct form is given here where you see that the "LoadDerivativeCoefficients" has an opposite sign than what is given in the InitializePDECoefficients ref page. So you changing the sign is the right thing to do. Note that $\alpha$ also has the wrong sign, but that is irrelevant in your example. I fixed the documentation a few weeks back when I realized the typo.

Concerning the boundary condition, note that we solve for the increment dU (not U) The boundary conditions, however, are made for applying to U. So we proceed by setting up the initial guess in such a way that it has the boundary values set, and at the same time modify the stiffness matrix to respect that in every iteration. So the increment of the boundary condition is 0.

| improve this answer | |
  • $\begingroup$ Thank you for pointing that out. I update my post based on your answer. I'm not sure that I understand it correctly or not. Please suggest. Thank you again.@user21 $\endgroup$ – Wilhelm Oct 24 '17 at 5:12
  • $\begingroup$ @Wilhelm, yes, I think you got that right. $\endgroup$ – user21 Oct 24 '17 at 13:20
  • $\begingroup$ Could you please remind me where can I find the definition of ProcessPDEEquations? I didn't find it in the Documentation. $\endgroup$ – Wilhelm Oct 26 '17 at 21:18

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.