Solving 1+u''(x) = b(x) with DiffEqOperators

Suppose b(x) is a function of a variable x, and let Dxx represent the second derivative d^2x/dx^2. Then the differential equation
1+Dxx u = b
has solution
u = (1+Dxx)\b.
I can’t figure out how to do this with DiffEqOperators.jl.

Let’s first consider the simpler equation: Dxx u = b. This can be solved as:

using DiffEqOperators

L = 2*pi
nx = 100
dx = L/(nx+1)

xpoints = range(dx, step=dx, length=nx)
b = sin.(xpoints)

Dxx = CenteredDifference(2, 2, dx, nx)
bc = Dirichlet0BC(Float64)

sol = (Dxx*bc)\b

Based on this I want to do something like (1+Dzz*bc)\b. However I can’t see how to make this work. I’ve tried using I from the LinearAlgebra package, the identity matrix, and a CenteredDifference of order zero, but all gave an error. Does anyone know what the solution is?

1 Like

Based on this question, I think I got it to work via:

using SparseArrays, LinearAlgebra
Axx,bxx = sparse(Dxx*bc)
Dh = I+Axx

and then doing Dh\b. Is this the best way to go about it?

The derivative operators are internally converted to sparse matrices anyway, so your approach is in principle good.

Beware that one needs to solve (I + Dxx*bc) \ (b-bxx) though. See for a discussion of ghost nodes and boundary conditions, and take a look at how DiffEqOperators implements \.

1 Like