🟢 Temperature model

This is a classical example from the Trilinos library.

This is a simple example in which we aim at solving $\Delta T+\alpha N(T,\beta)=0$ with boundary conditions $T(0) = T(1)=\beta$. This example is coded in examples/chan.jl. We start with some imports:

using BifurcationKit, Plots
const BK = BifurcationKit

N(x; a = 0.5, b = 0.01) = 1 + (x + a*x^2)/(1 + b*x^2)

We then write our functional:

function F_chan(x, p)
	(;α, β) = p
	f = similar(x)
	n = length(x)
	f[1] = x[1] - β
	f[n] = x[n] - β
	for i=2:n-1
		f[i] = (x[i-1] - 2 * x[i] + x[i+1]) * (n-1)^2 + α * N(x[i], b = β)
	end
	return f
end

We want to call a Newton solver. We first need an initial guess:

n = 101
sol0 = [(i-1)*(n-i)/n^2+0.1 for i=1:n]

# set of parameters
par = (α = 3.3, β = 0.01)

Finally, we need to provide some parameters for the Newton iterations. This is done by calling

optnewton = NewtonPar(tol = 1e-9, max_iterations = 10)

We call the Newton solver:

prob = BifurcationProblem(F_chan, sol0, par, (@optic _.α),
	# function to plot the solution
	plot_solution = (x, p; k...) -> plot!(x; ylabel="solution", label="", k...))
# we set verbose to true to see the newton iterations
sol = @time BK.solve(prob, Newton(), @set optnewton.verbose = true)

┌─────────────────────────────────────────────────────┐
│ Newton step         residual      linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│       0     │       2.3440e+01     │        0       │
│       1     │       1.3774e+00     │        1       │
│       2     │       1.6267e-02     │        1       │
│       3     │       2.4336e-06     │        1       │
│       4     │       6.7452e-12     │        1       │
└─────────────┴──────────────────────┴────────────────┘
  0.001050 seconds (384 allocations: 2.299 MiB)

Note that, in this case, we did not give the Jacobian. It was computed internally using Automatic Differentiation.

We can perform numerical continuation w.r.t. the parameter $\alpha$. This time, we need to provide additional parameters, but now for the continuation method:

optcont = ContinuationPar(max_steps = 150,
	p_min = 0., p_max = 4.2,
	newton_options =optnewton)

Next, we call the continuation routine as follows.

br = continuation(prob, PALC(), optcont; plot = true)
nothing #hide

The parameter axis lens = @optic _.α is used to extract the component of par corresponding to α. Internally, it is used as get(par, lens) which returns 3.3.

You should see

scene = title!("") #hide

The left figure is the norm of the solution as function of the parameter $p=\alpha$, the y-axis can be changed by passing a different record_from_solution to BifurcationProblem. The top right figure is the value of $\alpha$ as function of the iteration number. The bottom right is the solution for the current value of the parameter. This last plot can be modified by changing the argument plot_solution to BifurcationProblem.

What if we want to compute to continue both ways in one call?

br = continuation(prob, PALC(), optcont; bothside = true)
plot(br)