Bratu problem with contique

[adtzlr]

Hi @nschloe, @gdmcbain,

I just wanted to show you that I did something similar - though not as mature - to pacopy, called contique. It uses a magnitude of component - based numeric continuation of equilibrium curves (think of it as an extension to your natural method but not sticking to lambda in continuation).

I also tried to solve the Bratu problem and it looks good I think:

It automatically approximates the jacobian by finite differences if no jacobian of the problem is passed. So the Bratu problem is just:

def fun(x, lpf):
    n = len(x)
    h = 1 / (n - 1)
    A = np.diag(2 * np.ones_like(x) / h ** 2)
    for i in [1, -1]:
        A -= np.diag(np.ones_like(x[:-1]) / h ** 2, i)
    f = -A.dot(x) + lpf * np.exp(x)
    for i in [0, -1]:
        f[i] = x[i]
    return f

along with the definition of H and the calculation of the unknown u:

def val(x):
    n = len(x)
    h = 1 / (n - 1)
    H = np.ones(n) * h
    for i in [0, -1]:
        H[i] = h / 2
    return np.sqrt(x.dot(H * x))

...and passed to contique.solve:

# initial solution
n = 51
lpf0 = 0.0
x0 = np.zeros(n)

# numeric continuation
Res = contique.solve(
	fun=fun,
	x0=x0,
	lpf0=lpf0,
	dxmax=0.5,
	dlpfmax=0.5,
	maxsteps=22,
	tol=1e-10,
	rebalance=True,
	minlastfailed=0,
)

Results may be extracted with:

# extract/calculate load-proportionality factor and unknown u
lpf = np.array([res.x[-1] for res in Res])
u = np.array([val(res.x[:-1]) for res in Res])

# create figure
plt.plot(lpf, u, "x-")
plt.xlabel("load-proportionality-factor $LPF$")
plt.ylabel("$||u||_2$")
plt.plot(0, 0, "C0o", ms=10)
plt.savefig("test_bratu.svg")

Source-Code of the example is here: https://github.com/adtzlr/contique/blob/main/test/test_bratu.py

I thought of using pacopy beforce creating my own little module. But I wanted to code my own knowledge. Probably some pieces of my code could be interesting for you 😃

BR Andreas