{
"cells": [
{
"cell_type": "code",
"source": [
"using Revise\n",
"using PseudoArcLengthContinuation, LinearAlgebra, Plots, SparseArrays;\n",
"const Cont = PseudoArcLengthContinuation\n"
],
"outputs": [
{
"output_type": "execute_result",
"execution_count": 14,
"data": {
"text/plain": [
"PseudoArcLengthContinuation"
]
},
"metadata": {}
}
],
"execution_count": 14,
"metadata": {
"collapsed": false,
"outputHidden": false,
"inputHidden": false
}
},
{
"cell_type": "code",
"source": [
"f1(u, v) = u^2*v\n",
"\n",
"function F_bru(x, α, β; D1 = 0.008, D2 = 0.004, l = 1.0)\n",
"\tn = div(length(x), 2)\n",
"\th = 1.0 / (n+1); h2 = h*h\n",
"\n",
"\tu = @view x[1:n]\n",
"\tv = @view x[n+1:2n]\n",
"\n",
"\t# output\n",
"\tf = similar(x)\n",
"\n",
"\tf[1] = u[1] - α\n",
"\tf[n] = u[n] - α\n",
"\tfor i=2:n-1\n",
"\t\tf[i] = D1/l^2 * (u[i-1] - 2u[i] + u[i+1]) / h2 - (β + 1) * u[i] + α + f1(u[i], v[i])\n",
"\tend\n",
"\n\n",
"\tf[n+1] = v[1] - β / α\n",
"\tf[end] = v[n] - β / α;\n",
"\tfor i=2:n-1\n",
"\t\tf[n+i] = D2/l^2 * (v[i-1] - 2v[i] + v[i+1]) / h2 + β * u[i] - f1(u[i], v[i])\n",
"\tend\n",
"\n",
"\treturn f\n",
"end\n",
"\n",
"function Jac_mat(x, α, β; D1 = 0.008, D2 = 0.004, l = 1.0)\n",
"\tn = div(length(x), 2)\n",
"\th = 1.0 / (n+1); hh = h*h\n",
"\n",
"\tJ = zeros(2n, 2n)\n",
"\n",
"\tJ[1, 1] = 1.0\n",
"\tfor i=2:n-1\n",
"\t\tJ[i, i-1] = D1 / hh/l^2\n",
"\t\tJ[i, i] = -2D1 / hh/l^2 - (β + 1) + 2x[i] * x[i+n]\n",
"\t\tJ[i, i+1] = D1 / hh/l^2\n",
"\t\tJ[i, i+n] = x[i]^2\n",
"\tend\n",
"\tJ[n, n] = 1.0\n",
"\n",
"\tJ[n+1, n+1] = 1.0\n",
"\tfor i=n+2:2n-1\n",
"\t\tJ[i, i-n] = β - 2x[i-n] * x[i]\n",
"\t\tJ[i, i-1] = D2 / hh/l^2\n",
"\t\tJ[i, i] = -2D2 / hh/l^2 - x[i-n]^2\n",
"\t\tJ[i, i+1] = D2 / hh/l^2\n",
"\tend\n",
"\tJ[2n, 2n] = 1.0\n",
"\treturn J\n",
"end\n",
"\n",
"function Jac_sp(x, α, β; D1 = 0.008, D2 = 0.004, l = 1.0)\n",
"\t# compute the Jacobian using a sparse representation\n",
"\tn = div(length(x), 2)\n",
"\th = 1.0 / (n+1); hh = h*h\n",
"\n",
"\tdiag = zeros(2n)\n",
"\tdiagp1 = zeros(2n-1)\n",
"\tdiagm1 = zeros(2n-1)\n",
"\n",
"\tdiagpn = zeros(n)\n",
"\tdiagmn = zeros(n)\n",
"\n",
"\tdiag[1] = 1.0\n",
"\tdiag[n] = 1.0\n",
"\tdiag[n + 1] = 1.0\n",
"\tdiag[end] = 1.0\n",
"\n",
"\tfor i=2:n-1\n",
"\t\tdiagm1[i-1] = D1 / hh/l^2\n",
"\t\tdiag[i] = -2D1 / hh/l^2 - (β + 1) + 2x[i] * x[i+n]\n",
"\t\tdiagp1[i] = D1 / hh/l^2\n",
"\t\tdiagpn[i] = x[i]^2\n",
"\tend\n",
"\n",
"\tfor i=n+2:2n-1\n",
"\t\tdiagmn[i-n] = β - 2x[i-n] * x[i]\n",
"\t\tdiagm1[i-1] = D2 / hh/l^2\n",
"\t\tdiag[i] = -2D2 / hh/l^2 - x[i-n]^2\n",
"\t\tdiagp1[i] = D2 / hh/l^2\n",
"\tend\n",
"\treturn spdiagm(0 => diag, 1 => diagp1, -1 => diagm1, n => diagpn, -n => diagmn)\n",
"end\n",
"\n\n\n",
"function finalise_solution(z, tau, step, contResult)\n",
"\tn = div(length(z), 2)\n",
"\tprintstyled(color=:red, \"--> Solution constante = \", norm(diff(z[1:n])), \" - \", norm(diff(z[n+1:2n])), \"\\n\")\n",
"end"
],
"outputs": [
{
"output_type": "execute_result",
"execution_count": 15,
"data": {
"text/plain": [
"finalise_solution (generic function with 1 method)"
]
},
"metadata": {}
}
],
"execution_count": 15,
"metadata": {
"collapsed": false,
"outputHidden": false,
"inputHidden": false
}
},
{
"cell_type": "code",
"source": [
"n = 101\n",
"# const Δ = spdiagm(0=>2ones(N), -1=>0ones(N-1), 1=>-ones(N-1))\n",
"Jac_fd(u0, α, β, l = l) = Cont.finiteDifferences(u->F_bru(u, α, β, l=l), u0)\n",
"\n",
"a = 2.\n",
"b = 5.45\n",
"sol0 = vcat(a*ones(n), b/a*ones(n))\n",
"\n",
"opt_newton = Cont.NewtonPar(tol = 1e-11, verbose = true)\n",
"\t# ca fait dans les 60.2k Allocations\n",
"\tout, hist, flag = @time Cont.newton(\n",
"\t\tx -> F_bru(x, a, b),\n",
"\t\tx -> Jac_mat(x, a, b),\n",
"\t\tsol0,\n",
"\t\topt_newton)"
],
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"\n",
" Newton Iterations \n",
" Iterations Func-count f(x) Linear-Iterations\n",
"\n",
" 0 1 0.0000e+00 0\n",
" 0.137356 seconds (406.39 k allocations: 19.447 MiB, 7.24% gc time)\n"
]
},
{
"output_type": "execute_result",
"execution_count": 16,
"data": {
"text/plain": [
"([2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0 … 2.725, 2.725, 2.725, 2.725, 2.725, 2.725, 2.725, 2.725, 2.725, 2.725], [0.0], true, 0)"
]
},
"metadata": {}
}
],
"execution_count": 16,
"metadata": {
"collapsed": false,
"outputHidden": false,
"inputHidden": false
}
},
{
"cell_type": "code",
"source": [
"opts_br0 = ContinuationPar(dsmin = 0.001, dsmax = 0.0061, ds= 0.0051, pMax = 1.2, save = false, theta = 0.01, detect_fold = true, detect_bifurcation = true, nev = 16, plot_every_n_steps = 50)\n",
"\topts_br0.newtonOptions.maxIter = 20\n",
"\topts_br0.newtonOptions.tol = 1e-8\n",
"\topts_br0.maxSteps = 130\n",
"\n",
"\tbr, u1 = @time Cont.continuation(\n",
"\t\t(x, p) -> F_bru(x, a, b, l = p),\n",
"\t\t(x, p) -> Jac_mat(x, a, b, l = p),\n",
"\t\tout,\n",
"\t\t0.3,\n",
"\t\topts_br0,\n",
"\t\tplot = true,\n",
"\t\tplotsolution = (x;kwargs...)->(N = div(length(x), 2);plot!(x[1:N], subplot=4, label=\"\");plot!(x[N+1:2N], subplot=4, label=\"\")),\n",
"\t\tfinaliseSolution = finalise_solution,\n",
"\t\tprintsolution = x->norm(x, Inf64), verbosity = 0)\n"
],
"outputs": [
{
"output_type": "display_data",
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {}
},
{
"output_type": "stream",
"name": "stdout",
"text": [
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n",
"\u001b[31m--> Solution constante = 0.0 - 0.0\u001b[39m\n"
]
},
{
"output_type": "display_data",
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"image/svg+xml": [
"\n",
"