|
| 1 | +import numpy as np |
| 2 | +import scipy.optimize as opt |
| 3 | +import matplotlib.pylab as plt |
| 4 | + |
| 5 | +from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right |
| 6 | + |
| 7 | + |
| 8 | +def main(): |
| 9 | + def rho(x): |
| 10 | + return max(abs(np.linalg.eigvals(np.eye(M) - np.diag([x[i] for i in range(M)]).dot(coll.Qmat[1:, 1:])))) |
| 11 | + |
| 12 | + M = 2 |
| 13 | + |
| 14 | + coll = CollGaussRadau_Right(M, 0, 1) |
| 15 | + |
| 16 | + x0 = np.ones(M) |
| 17 | + d = opt.minimize(rho, x0, method='Nelder-Mead') |
| 18 | + |
| 19 | + print(d) |
| 20 | + numsteps = 800 |
| 21 | + xdim = np.linspace(0, 8, numsteps) |
| 22 | + ydim = np.linspace(0, 13, numsteps) |
| 23 | + |
| 24 | + minfield = np.zeros((len(xdim), len(ydim))) |
| 25 | + |
| 26 | + for idx, x in enumerate(xdim): |
| 27 | + for idy, y in enumerate(ydim): |
| 28 | + minfield[idx, idy] = max(abs(np.linalg.eigvals(np.eye(M) - np.diag([x, y]).dot(coll.Qmat[1:, 1:])))) |
| 29 | + |
| 30 | + # Set up plotting parameters |
| 31 | + params = {'legend.fontsize': 20, |
| 32 | + 'figure.figsize': (12, 8), |
| 33 | + 'axes.labelsize': 20, |
| 34 | + 'axes.titlesize': 20, |
| 35 | + 'xtick.labelsize': 16, |
| 36 | + 'ytick.labelsize': 16, |
| 37 | + 'lines.linewidth': 3 |
| 38 | + } |
| 39 | + plt.rcParams.update(params) |
| 40 | + |
| 41 | + plt.figure() |
| 42 | + plt.pcolor(xdim, ydim, minfield.T, cmap='Reds', vmin=0, vmax=1) |
| 43 | + plt.text(d.x[0], d.x[1], 'X', horizontalalignment='center', verticalalignment='center') |
| 44 | + plt.xlim((min(xdim), max(xdim))) |
| 45 | + plt.ylim((min(ydim), max(ydim))) |
| 46 | + plt.xlabel('component 1') |
| 47 | + plt.ylabel('component 2') |
| 48 | + cbar = plt.colorbar() |
| 49 | + cbar.set_label('spectral radius') |
| 50 | + |
| 51 | + fname = 'data/parallelSDC_minimizer_full.png' |
| 52 | + plt.savefig(fname, rasterized=True, bbox_inches='tight') |
| 53 | + |
| 54 | + plt.figure() |
| 55 | + xdim_part = xdim[int(0.25*numsteps):int(0.75*numsteps)+1] |
| 56 | + ydim_part = ydim[0:int(0.25*numsteps)] |
| 57 | + minfield_part = minfield[int(0.25*numsteps):int(0.75*numsteps)+1, 0:int(0.25*numsteps)] |
| 58 | + plt.pcolor(xdim_part, ydim_part, minfield_part.T, cmap='Reds', vmin=0, vmax=1) |
| 59 | + plt.text(d.x[0], d.x[1], 'X', horizontalalignment='center', verticalalignment='center') |
| 60 | + plt.xlim((min(xdim_part), max(xdim_part))) |
| 61 | + plt.ylim((min(ydim_part), max(ydim_part))) |
| 62 | + plt.xlabel('component 1') |
| 63 | + plt.ylabel('component 2') |
| 64 | + cbar = plt.colorbar() |
| 65 | + cbar.set_label('spectral radius') |
| 66 | + |
| 67 | + fname = 'data/parallelSDC_minimizer_zoom.png' |
| 68 | + plt.savefig(fname, rasterized=True, bbox_inches='tight') |
| 69 | + |
| 70 | + |
| 71 | +if __name__ == "__main__": |
| 72 | + main() |
0 commit comments