ApproxFun with Robin boundary conditions

New to Julia
1

Hi,

I am new to Julia and am interested in using the ApproxFun package to solve 1D boundary value problems. I found interesting examples here : ApproxFunExamples/PDEs/Manipulate ODEs and PDEs.ipynb at master · JuliaApproximation/ApproxFunExamples · GitHub but couldn’t find examples with Robin boundary conditions for instance. Let’s say I want to solve the following:

-u’‘+u=0
u’(0)=0, u’(1)+u(1)=1

How to implement that using ApproxFun?

Many thanks in advance.

1 Like
2

Hi,

I know this doesn’t directly answer your question, but BSplineKit.jl allows Robin-type boundary conditions.

The first thing you would need to do is to convert your boundary conditions to homogeneous BCs, for instance via the change of variables v = u - 1. Your equation then becomes:

\begin{align*} -v'' + v &= -1 \\ v'(0) &= 0 \\ v'(1) + v(1) &= 0 \end{align*}

The interface for solving boundary value problems with BSplineKit is currently not exactly high level. I can explain later how this works, but right now one would do it as follows:

using BSplineKit
using LinearAlgebra

ts = range(0, 1; length = 101)
B = BSplineBasis(BSplineOrder(4), ts)

# Homogeneous boundary conditions
bc_left = Derivative(1)
bc_right = Derivative(0) + Derivative(1)  # Robin-type BC

# Create recombined B-spline basis, such that each element of the basis individually satisfies
# the specified homogeneous BCs
R = RecombinedBSplineBasis(B, bc_left, bc_right)

M″ = galerkin_matrix(R, (Derivative(0), Derivative(2)))
M = galerkin_matrix(R, (Derivative(0), Derivative(0)))

A = M - M″
y = galerkin_projection(x -> -one(x), R)  # project right-hand side

vcoefs = A \ y  # solve linear system for basis coefficients

v = Spline(R, vcoefs)  # finally, create spline with the solution (can be evaluated, plotted, differentiated, ...)

## Verification

v′ = Derivative() * v
v″ = Derivative() * v′

v′(0)  # = 0
v(1) + v′(1)  # ≈ 0

# Verify graphically that v″ = v + 1

using CairoMakie

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = L"x")
lines!(ax, 0..1, x -> v(x) + 1; label = L"v(x) + 1")
lines!(ax, 0..1, x -> v″(x); linestyle = :dash, label = L"v″(x)")
axislegend(ax; position = (0, 1), labelsize = 20)
fig

solution
solution2400×1800 90.9 KB

1 Like
3

There is an example of solving a (nonlinear) boundary value problem with mixed boundary conditions here:

4

Have you tried a sum of a Dirichlet and a Neumann operator, or a sum of Evaluation operators? I think this should work.

5

Hi,

Thank you everyone for your answers.

@jipolanco: I will see if I can use that for my actual problem, although I would like to use ApproxFun in priority.

To Sevi: My problem is a linear one, thus I don’t want to use Newton, but I will give it a try.

@jishnub: Actually, I tried something like this, but It didn’t work (I am new to Julia so I don’t know the intricacies of these operators yet). I tried for instance

d = Domain(0…1)
L = -D^2-I
u = [Neumann(d)+(x>0.5)*Dirichlet(d);L] \ [[0,1],0]

but that didn’t give the right solution. How would you impose Robin boundary conditions using a sum of Dirichlet and Neumann operators?

Thanks in advance.

6 
julia> d = Domain(0..1);

julia> D = Derivative(Chebyshev(d));

julia> L = -D^2-I;

julia> f = [ldirichlet(); rdirichlet() + rneumann(); L] \ [0, 1, 0];

julia> f(0)
0.0

julia> f(1) + f'(1)
1.0
7

Thank you very much @jishnub. Is there an option to add somewhere as I’m interested in a high accuracy approximation? Thanks in advance.

8

The operator \ accepts a tolerance, so instead of A \ b, you may call it as (A, b, tolerance=1e-16) (or whatever tolerance you want).

If you want higher precision, perhaps you could try using a BigFloat or a BigInt domain? This should also automatically bump the tolerances used. Something like d = Domain(big(0)..big(1)).