🟢 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 Plotsimport BifurcationKit as BKimport BifurcationKit: @optic, @setN(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 fend

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

n = 101sol0 = [(i-1)*(n-i)/n^2+0.1 for i=1:n]# set of parameterspar = (α = 3.3, β = 0.01)

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

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

We call the Newton solver:

prob = BK.ODEBifProblem(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 iterationssol = @time BK.solve(prob, BK.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.000916 seconds (317 allocations: 2.296 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 = BK.ContinuationPar(max_steps = 150,	p_min = 0., p_max = 4.2,	newton_options = optnewton)

Next, we call the continuation routine as follows.

br = BK.continuation(prob, BK.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

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 = BK.continuation(prob, BK.PALC(), optcont; bothside = true)plot(br)