An example of a boundary-value problem.

In a boundary value problem, we look for a function that satisfies a certain property.

Here, we look at solving the Boundary Value Problem: $$ u''(x) = f(x) $$ where x ranges between $[0, 1]$ and $u(0) = \alpha$ and $u(1) = \beta$.

This is different from a standard ODE because we have TWO constraints, one at $u(0)$ and one at $u(1)$.

To solve use a finite difference discretization on a grid.

We use a grid presentation of the function $u(x)$ on the interval $[0,1]$ so we divide $[0,1]$ up into $N+1$ increments of $1/N$ each.

This means that u is a vector with n = N+1 elements, but two of them are fixed at alpha and beta.

In [2]:
using Plots
N = 10
xgrid = collect(0:1/N:1)

n = length(xgrid)

example_u = zeros(n);
example_u[1] = 1       # suppose alpha = 1
example_u[end] = 1/2   # suppose beta = 1/2

example_u[2:end-1] = 1.06-1*(xgrid[2:end-1]-0.25).^2;

plot(xgrid, example_u)
[Plots.jl] Initializing backend: pyplot
WARNING: No working GUI backend found for matplotlib.
Out[2]:
In [5]:
N = 10
xgrid = collect(0:1/N:1)

alpha = 1
beta = 0.5
f = x -> -2.0
h = 1/N

A = zeros(N+1,N+1)
b = zeros(N+1,1)
for i=1:N+1
    if i == 1
        A[i,i] = 1
        b[i] = alpha
    elseif i==N+1
        A[i,i] = 1
        b[i] = beta
    else
        A[i,i-1] = 1/h^2
        A[i,i] = (-2)/h^2
        A[i,i+1] = 1/h^2
        b[i] = f(xgrid[i])
    end
end
u = A\b
plot(xgrid, u, marker=:circle)
Out[5]:
In [31]:
plot(xgrid, u, legend=false, marker=:circle, linewidth=1.25, markersize=6,size=(100*3,100*3))
xlabel!("x")
ylabel!("u")
xticks!(0:.1:1)
yticks!(0:.1:1.2)
ylims!(0,1.25)
PyPlot.gcf()[:set_size_inches](3,3)
savefig("bvp-example.pdf")
In [6]:
for i=1:N+1
    for j=1:N+1
        @printf("%i & ", round(A[i,j]))
    end
    @printf(" \\\\ \n")
end
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &  \\ 
100 & -200 & 100 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &  \\ 
0 & 100 & -200 & 100 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &  \\ 
0 & 0 & 100 & -200 & 100 & 0 & 0 & 0 & 0 & 0 & 0 &  \\ 
0 & 0 & 0 & 100 & -200 & 100 & 0 & 0 & 0 & 0 & 0 &  \\ 
0 & 0 & 0 & 0 & 100 & -200 & 100 & 0 & 0 & 0 & 0 &  \\ 
0 & 0 & 0 & 0 & 0 & 100 & -200 & 100 & 0 & 0 & 0 &  \\ 
0 & 0 & 0 & 0 & 0 & 0 & 100 & -200 & 100 & 0 & 0 &  \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 100 & -200 & 100 & 0 &  \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 100 & -200 & 100 &  \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 &  \\ 

An example of the steady-state heat equation

Here, we look at solving the Boundary Value Problem:

$$ (\partial ( p(x) ( \partial u/\partial x ) ) / \partial x = f(x) $$

where x ranges between [0, 1] and u(0) = alpha and u(1) = beta.

This means we are looking for a function, u(x) where

$$(p(x) u'(x))' = f(x)$$$$ p(x) u''(x) + p'(x) u'(x) = f(x) $$

and u(0) = alpha and u(1) = beta.

Approximate the derivatives $$ (\partial ( p(x) ( \partial u/\partial x ) ) / \partial x = f(x) $$ outside-in. To make this clear, lets call $$ g(x) = p(x) ( \partial u/\partial x ) $$ In which case, we have partial $g(x) / \partial x = f(x)$ so if we use a central difference formula with width "h", we get $$ f(x) = 1/h ( g(x+h/2) - g(x-h/2) ) $$

But, then, we can use the same idea to evaluate g(x+h/2)

$$ g(x+h/2) = p(x+h/2) 1/h (u(x+h) - u(x)) $$$$ g(x-h/2) = p(x-h/2) 1/h (u(x) - u(x-h)) $$

Thus, we have a system of linear equations we can solve for f(x), and u(x). Suppose that $u_i = u(x_i), u_{i+1} = u(x_i + h), $u_{i-1} = u(x_i - h)$, then:

In [39]:
# The solution strategy
# The idea with BVPs is to approximate the solution on a grid of points.
# Let's divide the region into 11 values

N = 10
xgrid = collect(0:1/N:1)

alpha = -1
beta = 1
f = x -> -2*sin(pi*x)   # we are heating a rod in the middle
p = x -> 1-4*(x-0.5)^2  # the rod becomes more like an insulator towards the end.
h = 1/N

A = zeros(N+1,N+1)
b = zeros(N+1,1)
for i=1:N+1
    if i == 1
        A[i,i] = 1
        b[i] = alpha
    elseif i==N+1
        A[i,i] = 1
        b[i] = beta
    else
        A[i,i-1] = p(xgrid[i] - h/2)/h^2
        A[i,i] = (-p(xgrid[i] + h/2)-p(xgrid[i] - h/2))/h^2
        A[i,i+1] = p(xgrid[i] + h/2)/h^2
        b[i] = f(xgrid[i])
    end
end
u = A\b
plot(xgrid, u, marker=:circle)
Out[39]: