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

General Usage
1

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
2

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?

3

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 https://github.com/JuliaDiffEq/DiffEqOperators.jl/files/3267835/ghost_node.pdf for a discussion of ghost nodes and boundary conditions, and take a look at how DiffEqOperators implements \.

1 Like