๐ข 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, LinearAlgebra, Plots, Parameters
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)
@unpack ฮฑ, ฮฒ = 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, (@lens _.ฮฑ),
# 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 newton( prob, @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.001113 seconds (344 allocations: 1.262 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 = @lens _.ฮฑ
is used to extract the component of par
corresponding to ฮฑ
. Internally, it is used as get(par, lens)
which returns 3.3
.
We don't need to call newton
first in order to use continuation
.
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 recordFromSolution
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 plotSolution
to BifurcationProblem
.
Two Fold points were detected. This can be seen by looking at br.specialpoint
, by the black dots on the continuation plots when doing plot(br, plotfold=true)
or by typing br
in the REPL. Note that the bifurcation points are located in br.specialpoint
.
What if we want to compute to continue both ways in one call?
br = continuation(prob, PALC(), optcont; bothside = true)
plot(br)