restart; read `Lanczos.map`; a := 0: b := infinity: omega := t -> exp (-t): iprod := (f, g, x) -> int (f * g * omega (x), x=a..b): n := 10: Lanczos (iprod, 'alpha', 'beta', n): print (alpha, beta); read `Gram.map`: read `SackDonovan.map`: N := 10: phi := 1: alphahat := array (1..2*N-1, [seq (2*k-1, k=1..2*N-1)]): betahat := array (1..2*N-1, [seq (k, k=1..2*N-1)]): gammahat := array (1..2*N-2, [seq (k, k=1..2*N-2)]): L := (n, x) -> sum ('(-1)^(n-k) * binomial (n,k) * x^k / k!', 'k'=0..n): nucos[0] := 1 + 1 / (1+phi^2): for n from 1 to 2*N-1 do nucos[n] := (-1)^n * phi^n / (1+phi^2)^(n+1) * sum ('binomial (n+1,2*k+1) * (-1)^k * phi^(n-2*k)', 'k'=0..floor (n/2)); od: nusin[0] := 1 + phi / (1+phi^2): for n from 1 to 2*N-1 do nusin[n] := (-1)^n * phi^n / (1+phi^2)^(n+1) * sum ('binomial (n+1,2*k) * (-1)^k * phi^(n+1-2*k)', 'k'=0..floor ((n+1)/2)); od: SackDonovan (L (0, x), alphahat, betahat, gammahat, nucos, N, 'alphacos', 'betacos'); SackDonovan (L (0, x), alphahat, betahat, gammahat, nusin, N, 'alphasin', 'betasin'); print (alphacos, betacos); print (alphasin, betasin); read `Lobatto.map`: N := 11: Digits := 100: a := 1: b := 2: omega := x -> 1 / x: alphahat := array (1..2*N-1, [0$2*N-1]): betahat := array (1..2*N-1, [1$2*N-1]): gammahat := array (1..2*N-2, [0$2*N-2]): pi0 := 1: nu[0] := ln (b) - ln (a): for n from 1 to 2*N-1 do nu[n] := (b^n - a^n) / n; od: mu0 := nu[0] / pi0: SackDonovan (pi0, alphahat, betahat, gammahat, nu, N, 'alpha', 'beta'); Lobatto (mu0, evalf(op(alpha)), evalf(op(beta)), N-1, a, b, 'x', 'w'); print (x, w);