# Solving non-linear partial differential equation numerically: $u_{xx}+u_{yy}=\mathrm{e}^{u}$

To start with, I need to solve this partial equation numerically, but I do not know how to do that. If I try a finite difference method, I face a problem that $$u_{i,j}$$ is also located in exponential, consequently, there is no way to find $$u_{i,j}$$. I tried this fragmentation: $$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta y^2} = \mathrm{e}^{u_{i,j}}$$ What can help me to solve this equation? What ways do exist to solve this problems?

Boundary condiditions: $$u(0,y)=a, \quad u_x(0,y) = b$$ $$u(x,0)=c, \quad u_y(x,0) = d$$

• For a beginner-level approach, one could add the time-derivative term there and solve it as a time-evolution problem to the steady state, using the method of lines and a standard ODE integrator. Commented Nov 22, 2023 at 23:34
• You just happen to get a nonlinear system of equations, as shown in your first formula. You solve that like any other nonlinear system, for example with Newton's method. Commented Nov 23, 2023 at 0:35
• If you want more insight into the exact methods, solvers, and preconditioners that are used on this problem, it is called the Bratu problem and is a decent amount of literature on its numerical solution. Commented Nov 23, 2023 at 23:06
• You are using two boundary conditions on $x=0$ and two conditions on $y=0$, and nothing on the other parts of the boundary. Is this what you want to do ? Commented Nov 28, 2023 at 13:28

You can use something like MethodOfLines.jl to solve your problem. You can modify the following example, something like this

using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve

@parameters x y
@variables u(..)

Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ exp(u(x,y))

a = 5.2
b = 4.3
c = 3.5
d = 2.2

bcs = [u(0, y) ~ a,
u(x, 0) ~ c,
Dx(u(0, y)) ~ b,
Dy(u(x, 0)) ~ d]

domains = [x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])

dx = 0.1
dy = 0.1

discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)

prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson())

u_sol = sol[u(x, y)]

using Plots

heatmap(sol[x], sol[y], u_sol, xlabel="x", ylabel="y",
title="Solution")


Where I've plugged in some random boundary conditions. You might want to tweak some parameters like tolerance, maximum iterations to better suit your problem.