I'm working on trying to solve a steady state PDE in python and I am quite new to python and I am slightly confused on the implementation details. I have an example implementation of solving a non-linear Laplace equation: $$-\nu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \ = \ f(u)$$ but I am unsure how I can edit this python file to solve for:

$$\left(\mu_{x}\frac{\partial^{2}u}{\partial x^{2}}+\mu_{y}\frac{\partial^{2}u}{\partial y^{2}}\right)=f(x,y) \tag 1$$

and you discretize the problem using second-order finite difference formulas, leading to the discretized form:

$$- \mu_{x} \left(\frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{h^2}\right) - \mu_{y} \left(\frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{h^2}\right) = f_{i,j}$$

I am trying to numerically solve the problem for the right-hand side by evaluating the function:

$$f(x, y) = \begin{cases} 1 & \text{if } x ≥ 0.1 and x ≤ 0.3 and y≥ 0.1 and y ≤ 0.3 \\ 0 & \text{otherwise} \\ \end{cases}$$

in python but I'm not entirely sure how, I have most of the basis for it down but I am unsure how to now implement this equation to begin solving.

I have provided my current code: solver.py which solves the nonlinear laplacian equation and I am wondering if someone could explain how I might edit this to solve my desired equation.

Boundary conditions clarifications: $$$$-(\mu_x \frac{\partial^2 u}{\partial x^2} + \mu_y \frac{\partial^2 u}{\partial y^2}) = f(x, y), \quad (y = 0, 1) \times (x, 0, 1),$$$$ $$$$\text{and the boundary condition } u = 0 \text{ on } \Gamma \text{ (the boundary of the unit square)} f(x,y) \text{ is a function modelling a heat source.} \mu_x > 0 \text{ and } \mu_y > 0 \text{ are heat conductivity heat coefficients},$$$$

• This question is not well written. In essence it reads "I implemented a start of something, but I don't know how to complete it." But you only link to a place where you show 13 lines of code, you don't say what your problem is, what you're stuck on, and what you already tried. Your post also does not have a question mark. Dec 9, 2023 at 0:13
• How did you derive the Laplacian matrix in your original code? Can you think about how you might modify that to apply to this problem? Dec 10, 2023 at 19:08
• @whpowell96 Thanks for the reply! I think the only thing that needs changing for the matrix section, when doing the new problem is the stencil? it should be: stencil = np.array([4., mu_x * -1., mu_x * -1., mu_y * -1., mu_y * -1.]) / h**2?
– blov
Dec 10, 2023 at 20:16
• I think the 4 needs to be something different that contains a combination of both coefficients. If you test this and it works, please write up an answer with the details so that future users may find this solution via the search function. Dec 11, 2023 at 0:29