I am attempting to follow this tutorial in the documentation on using FEM to solve PDEs. I am having difficulty understanding how to formulate the Neumann boundary condition for my free-boundary problem. In the notation of the tutorial, it seems integral to be able to write down the equation in the form (I drop the "time" terms since they are not relevant to my question) \begin{equation} \nabla \cdot (-c(t,X,u,\nabla_X u) \nabla u-\alpha (t,X,u,\nabla _Xu) u+\gamma (t,X,u,\nabla _Xu)) + \beta (t,X,u,\nabla _Xu)\cdot \nabla u+a(t,X,u,\nabla _Xu) u - f(t,X,u,\nabla _Xu) = 0\ , \end{equation}

since the term in the first set of parentheses is used to define the Neumann boundary condition:

$$ \overset{\rightharpoonup }{n}\cdot (c \nabla u-\gamma +\alpha u)=g-q u\ . $$

Well, my PDE for $F = F(p, w)$ is

$$ \frac{1}{2}\sigma^2 p^2 \partial_p^2 F + \mu p \partial p + (\rho w + \xi) \partial_w F + (p \partial_w F)^\gamma - \beta F = 0\ , $$

where all the lowercase Greek letters are constants.

(i) I am first of all unclear as to how best to re-write this in the standard form above, since I have a first-order derivative in $w$ but a second-order derivative in $p$. In fact, I can see multiple ways to re-write the equation, but am unsure which one is correct.

(ii) Once this has been accomplished, I would like the following "smooth-fit" conditions to hold on the boundary $\partial \Omega$ of the domain over which I solve the PDE:

$$ F = G,\ \quad \partial_p F = \partial_p G, \quad \partial_w F = \partial_w G. $$

Here $G(p, w)$ is a known function. The first condition is a Dirichlet condition which is straightforward to implement. How would I formulate the latter two conditions in the Mathematica syntax for the Neumann boundary condition?

  • $\begingroup$ I am having trouble displaying the equations: see this meta post: mathematica.meta.stackexchange.com/q/2746/7984 Edits welcome! $\endgroup$
    – Anthony
    Commented Jul 4, 2023 at 7:13
  • $\begingroup$ fixed you latex. do not add ``` and ``` around the math. $\endgroup$
    – Nasser
    Commented Jul 4, 2023 at 9:29
  • $\begingroup$ Have you tried solving this with the tensor grid product method? Might be easier then to convert this in standard form? $\endgroup$
    – user21
    Commented Jul 19, 2023 at 4:38
  • $\begingroup$ While it is probably possible to convert this into the standard coefficient form, it does not look like the Neumann values directly fit to that form. $\endgroup$
    – user21
    Commented Jul 19, 2023 at 4:40


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.