Coverage for pySDC/projects/AsympConv/conv_test_toinf.py: 0%

1import matplotlib 

2import numpy as np 

3import scipy.linalg as LA 

4 

5matplotlib.use('Agg') 

6import matplotlib.pylab as plt 

7from matplotlib import rc 

8 

9from pySDC.core.Collocation import CollBase 

10 

11 

12def compute_and_plot_specrad(Nnodes, lam): 

13 """ 

14 Compute and plot the spectral radius of the smoother for different step-sizes 

15 

16 Args: 

17 Nnodes: number of collocation nodes 

18 lam: test parameter representing the spatial problem 

19 """ 

20 

21 Nsteps = 1 

22 

23 coll = CollBase(Nnodes, 0, 1, quad_type='RADAU-RIGHT') 

24 Qmat = coll.Qmat[1:, 1:] 

25 

26 # do LU decomposition of QT (St. Martin's trick) 

27 QT = coll.Qmat[1:, 1:].T 

28 [_, _, U] = LA.lu(QT, overwrite_a=True) 

29 QDmat = U.T 

30 

31 Nmat = np.zeros((Nnodes, Nnodes)) 

32 Nmat[:, -1] = 1 

33 

34 Nsweep_list = [1, Nnodes - 1, Nnodes] 

35 color_list = ['red', 'blue', 'green'] 

36 marker_list = ['s', 'o', 'd'] 

37 

38 setup_list = zip(Nsweep_list, color_list, marker_list) 

39 

40 xlist = [10**i for i in range(11)] 

41 

42 rc('font', **{"sans-serif": ["Arial"], "size": 24}) 

43 plt.subplots(figsize=(15, 10)) 

44 

45 for Nsweeps, color, marker in setup_list: 

46 Emat = np.zeros((Nsteps, Nsteps)) 

47 np.fill_diagonal(Emat[1:, :], 1) 

48 

49 Prho_list = [] 

50 predict_list = [] 

51 for x in xlist: 

52 mat = np.linalg.inv(np.eye(Nnodes * Nsteps) - x * lam * np.kron(np.eye(Nsteps), QDmat)).dot( 

53 x * lam * np.kron(np.eye(Nsteps), (Qmat - QDmat)) + np.kron(Emat, Nmat) 

54 ) 

55 mat = np.linalg.matrix_power(mat, Nsweeps) 

56 

57 Prho_list.append(max(abs(np.linalg.eigvals(mat)))) 

58 predict_list.append(1.0 / x) 

59 

60 if len(predict_list) > 1: 

61 print( 

62 x, 

63 predict_list[-1], 

64 Prho_list[-1], 

65 Prho_list[-2] / Prho_list[-1], 

66 predict_list[-2] / predict_list[-1], 

67 ) 

68 

69 plt.loglog( 

70 xlist, 

71 Prho_list, 

72 linestyle='-', 

73 linewidth=3, 

74 color=color, 

75 marker=marker, 

76 markersize=10, 

77 label='spectral radius, L=' + str(Nsteps), 

78 ) 

79 

80 plt.loglog( 

81 xlist, 

82 [item / predict_list[0] for item in predict_list], 

83 linestyle='--', 

84 linewidth=2, 

85 color='k', 

86 label='estimate', 

87 ) 

88 

89 plt.xlabel('time-step size') 

90 plt.ylabel('spectral radius') 

91 plt.legend(loc=3, numpoints=1) 

92 plt.grid() 

93 plt.ylim([1e-16, 1e00]) 

94 

95 if type(lam) is complex: 

96 fname = 'data/smoother_specrad_toinf_M' + str(Nnodes) + '_LU_imag.png' 

97 else: 

98 fname = 'data/smoother_specrad_toinf_M' + str(Nnodes) + '_LU_real.png' 

99 plt.savefig(fname, transparent=True, bbox_inches='tight') 

100 

101 

102if __name__ == "__main__": 

103 compute_and_plot_specrad(Nnodes=3, lam=-1) 

104 compute_and_plot_specrad(Nnodes=3, lam=1j) 

105 compute_and_plot_specrad(Nnodes=7, lam=-1) 

106 compute_and_plot_specrad(Nnodes=7, lam=1j)