| 45 | |
| 46 | |
| 47 | def get_phi_psi(k, base): |
| 48 | |
| 49 | x = Symbol('x') |
| 50 | phi_coeff = np.zeros((k,k)) |
| 51 | phi_2x_coeff = np.zeros((k,k)) |
| 52 | if base == 'legendre': |
| 53 | for ki in range(k): |
| 54 | coeff_ = Poly(legendre(ki, 2*x-1), x).all_coeffs() |
| 55 | phi_coeff[ki,:ki+1] = np.flip(np.sqrt(2*ki+1) * np.array(coeff_).astype(np.float64)) |
| 56 | coeff_ = Poly(legendre(ki, 4*x-1), x).all_coeffs() |
| 57 | phi_2x_coeff[ki,:ki+1] = np.flip(np.sqrt(2) * np.sqrt(2*ki+1) * np.array(coeff_).astype(np.float64)) |
| 58 | |
| 59 | psi1_coeff = np.zeros((k, k)) |
| 60 | psi2_coeff = np.zeros((k, k)) |
| 61 | for ki in range(k): |
| 62 | psi1_coeff[ki,:] = phi_2x_coeff[ki,:] |
| 63 | for i in range(k): |
| 64 | a = phi_2x_coeff[ki,:ki+1] |
| 65 | b = phi_coeff[i, :i+1] |
| 66 | prod_ = np.convolve(a, b) |
| 67 | prod_[np.abs(prod_)<1e-8] = 0 |
| 68 | proj_ = (prod_ * 1/(np.arange(len(prod_))+1) * np.power(0.5, 1+np.arange(len(prod_)))).sum() |
| 69 | psi1_coeff[ki,:] -= proj_ * phi_coeff[i,:] |
| 70 | psi2_coeff[ki,:] -= proj_ * phi_coeff[i,:] |
| 71 | for j in range(ki): |
| 72 | a = phi_2x_coeff[ki,:ki+1] |
| 73 | b = psi1_coeff[j, :] |
| 74 | prod_ = np.convolve(a, b) |
| 75 | prod_[np.abs(prod_)<1e-8] = 0 |
| 76 | proj_ = (prod_ * 1/(np.arange(len(prod_))+1) * np.power(0.5, 1+np.arange(len(prod_)))).sum() |
| 77 | psi1_coeff[ki,:] -= proj_ * psi1_coeff[j,:] |
| 78 | psi2_coeff[ki,:] -= proj_ * psi2_coeff[j,:] |
| 79 | |
| 80 | a = psi1_coeff[ki,:] |
| 81 | prod_ = np.convolve(a, a) |
| 82 | prod_[np.abs(prod_)<1e-8] = 0 |
| 83 | norm1 = (prod_ * 1/(np.arange(len(prod_))+1) * np.power(0.5, 1+np.arange(len(prod_)))).sum() |
| 84 | |
| 85 | a = psi2_coeff[ki,:] |
| 86 | prod_ = np.convolve(a, a) |
| 87 | prod_[np.abs(prod_)<1e-8] = 0 |
| 88 | norm2 = (prod_ * 1/(np.arange(len(prod_))+1) * (1-np.power(0.5, 1+np.arange(len(prod_))))).sum() |
| 89 | norm_ = np.sqrt(norm1 + norm2) |
| 90 | psi1_coeff[ki,:] /= norm_ |
| 91 | psi2_coeff[ki,:] /= norm_ |
| 92 | psi1_coeff[np.abs(psi1_coeff)<1e-8] = 0 |
| 93 | psi2_coeff[np.abs(psi2_coeff)<1e-8] = 0 |
| 94 | |
| 95 | phi = [np.poly1d(np.flip(phi_coeff[i,:])) for i in range(k)] |
| 96 | psi1 = [np.poly1d(np.flip(psi1_coeff[i,:])) for i in range(k)] |
| 97 | psi2 = [np.poly1d(np.flip(psi2_coeff[i,:])) for i in range(k)] |
| 98 | |
| 99 | elif base == 'chebyshev': |
| 100 | for ki in range(k): |
| 101 | if ki == 0: |
| 102 | phi_coeff[ki,:ki+1] = np.sqrt(2/np.pi) |
| 103 | phi_2x_coeff[ki,:ki+1] = np.sqrt(2/np.pi) * np.sqrt(2) |
| 104 | else: |