CoCalc Public Filesmonic-sigmoidal.sagewsOpen with one click!
Author: Namig Guliyev
Views : 51
Compute Environment: Ubuntu 18.04 (Deprecated)
import sys def continued_fraction(r): a, b = floor(r), r - floor(r) value = [a] while b != 0: a, b = floor(1/b), 1/b - floor(1/b) value += [a] return value def positive_rational(n): """ Returns the n-th term of the Calkin-Wilf sequence """ b = bin(n)[2:] ld = '1' c = 0 cf = [] for i in range(len(b)-1, -1, -1): if b[i] == ld: c += 1 else: cf += [c] ld = b[i] c = 1 cf += [c] r = cf[-1] for k in range(len(cf)-2, -1, -1): r = cf[k] + 1 / r return r def rational_number(n): """ Returns the n-th rational number """ if n == 0: return 0 elif n % 2 == 0: return positive_rational(n // 2) else: return -rational_number(n+1) def polynomial(n): """ Returns the n-th monic polynomial as a list of coefficients """ if n == 1: return [1] else: cf = continued_fraction(positive_rational(n)) if len(cf) == 1: cf = [cf[0] - 2] else: cf = [cf[0]] + [k-1 for k in cf[1:-1]] + [cf[-1]-2] return [rational_number(k) for k in cf] + [1] def polynomial_function(p): """ list of coefficients --> function """ def u(t): value = 0 for i in range(len(p)): value += p[i] * (t ** i) return value return u beta1 = lambda t: exp(-1/t) if t > 0 else 0 beta = lambda a, b, t: beta1(b-t) / (beta1(b-t) + beta1(t-a)) # smooth transition h = lambda d, lambda_, t: 1 - min(0.5, lambda_) / (log(t - d + 1) + 1) if t > d else 1 - min(0.5, lambda_) coef_lin_map = lambda a, b, c, d: ((b*c - a*d) / (b - a), (d - c) / (b - a)) # coefficients of the linear map [a,b] -> [c,d] lower_bound = lambda p: p[0] + sum([c - abs(c) for c in p[1:]]) / 2 upper_bound = lambda p: p[0] + sum([c + abs(c) for c in p[1:]]) / 2 upper_bound_for_derivative = lambda p, a, b: sum([k * abs(p[k]) * (max(abs(a), abs(b)) ** (k-1)) for k in range(1, len(p))]) def sigma(d, lambda_): def sigma1(t): if t >= d: m = round(t / (2*d)) if t == (2*m+1)*d: # preventing possible rounding error m += 1 p = polynomial(m) um = polynomial_function(p) if t <= 2*m*d: M = h(d, lambda_, (2*m+1)*d) if len(p) == 1: # um is constant value = (1 + M) / 2 else: A1, A2 = lower_bound(p), upper_bound(p) am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3) value = am + bm * um(t / d - 2 * m + 1) else: K = (sigma1(2*m*d) + sigma1((2*m+1)*d)) / 2 if t - 2*m*d < d / 2: M = h(d, lambda_, (2*m+1)*d) if len(p) == 1: # um is constant v = (1 + M) / 2 delta = d / 2 else: A1, A2 = lower_bound(p), upper_bound(p) am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3) v = am + bm * um(t / d - 2 * m + 1) epsilon = (1 - M) / 6 C = upper_bound_for_derivative(p, 1, 1.5) delta = min(d / 2, (epsilon * d) / (bm * C)) value = K - (K - v) * beta(2*m*d, 2*m*d + delta, t) else: p = polynomial(m+1) um = polynomial_function(p) M = h(d, lambda_, (2*m+3)*d) if len(p) == 1: # um is constant v = (1 + M) / 2 delta = d / 2 else: A1, A2 = lower_bound(p), upper_bound(p) am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3) v = am + bm * um(t / d - 2 * m - 1) epsilon = (1 - M) / 6 C = upper_bound_for_derivative(p, -0.5, 0) delta = min(d / 2, (epsilon * d) / (bm * C)) value = K - (K - v) * (1 - beta((2*m+1)*d - delta, (2*m+1)*d, t)) else: value = (1 + h(d, lambda_, 3*d)) * (1 - beta1(d - t)) / 2 return value return sigma1 def anti_fusc(r): if r == 0: return 0 cf = continued_fraction(r) if len(cf) % 2 == 0: cf = cf[:-1] + [cf[-1]-1] + [1] value = 0 for i in range(len(cf)-1, -1, -1): try: value = value * (2 ** cf[i]) except: # Error: Position of the polynomial is too large to calculate return -2 if i % 2 == 0: value += (2 ** cf[i]) - 1 return value def anti_rational(r): """ Returns the index of r """ if r >= 0: return 2 * anti_fusc(r) else: return 2 * anti_fusc(-r) - 1 def anti_polynomial(p_): if p_ == [1]: return 1 del p_[-1] p = [anti_rational(r) for r in p_] if len(p) == 1: p = [p[0] + 2] else: p = [p[0]] + [i+1 for i in p[1:-1]] + [p[-1]+2] r = p[-1] for k in range(len(p)-2, -1, -1): try: r = p[k] + 1 / r except (MemoryError, OverflowError): # Error: Position of the polynomial is too large to calculate return 0 try: return anti_fusc(r) except (MemoryError, OverflowError): # Error: Position of the polynomial is too large to calculate return 0 def rational_approximating_number(x, epsilon): for k in range(1, ceil(1 / (2 * epsilon)) + 1): y = round(x * k) / k if abs(x - y) <= epsilon: return y def rational_approximating_Taylor_polynomial(f, a, b, epsilon): epsilon1 = epsilon / 2 n = 0 err = plot(abs(diff(f, x)), (x, a, b)).ymax() * (b - a) / 2 while (err >= epsilon1) & (n < 101): n += 1 err = plot(abs(diff(f, x, n+1)), (x, a, b)).ymax() * ((b - a) ** (n+1)) / (factorial(n+1) * (2 ** (n+1))) if (n == 101): return [0, 0] fp = taylor(f, x, (a + b) / 2, n).coefficients(sparse=False) tp = [] for k in range(len(fp)): tt = 0 for m in range(k, len(fp)): tt += fp[m] * (a ** (m-k)) * binomial(m, k) tp += [tt * ((b - a) ** k)] epsilon2 = epsilon - err p = [] for m in range(len(tp)-1, -1, -1): if m == 0: epsilon2 *= 2 r = rational_approximating_number(tp[m].n(), epsilon2 / 2) p = [r] + p epsilon2 -= abs(tp[m].n() - r.n()) while (p[-1] == 0) & (len(p) > 1): del p[-1] return p def coefficients_of_approximation(lambda_, f, a, b, epsilon): d = b - a p = rational_approximating_Taylor_polynomial(f, a, b, epsilon) if p == [0, 0]: return [0, 0, 0, 1] plc = p[-1] monic = [k / plc for k in p] + [1] if plc != 0 else [1] m = anti_polynomial(monic) if m <= 0: # Error! return [0, 0, 0, 0] M = h(d, lambda_, (2*m+1)*d) if m == 1: am, bm = 1 / 2, h(d, lambda_, 3*d) / 2 else: A1, A2 = lower_bound(monic), upper_bound(monic) am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3) c1, c2 = plc / bm, 2 * am * plc / (bm * (1 + h(d, lambda_, 3*d))) return [c1, c2, b - (2 * m * d), 2 * a - b] def graph_of_sigma(d, lambda_, xmin, xmax, text_ = ""): s = sigma(d, lambda_) h_ = lambda t: 1 - min(0.5, lambda_) / (log(t - d + 1) + 1) if t >= d else 1 - min(0.5, lambda_) image = plot(s, (x, xmin, xmax), color = 'blue') + plot(h_, (x, xmin, xmax), color = 'red', linestyle = "--") + plot(1, (x, xmin, xmax), color = 'green', linestyle = "--") if text_ == "lambda": image += text(r"$\lambda = " + str(lambda_).rstrip('0') + "$", (0.95 * xmax, 0.5 * image.ymin() + 0.5 * image.ymax()), fontsize = 15, rgbcolor = (0, 0, 0)) elif text_ == "s": image += text(r"$s = " + str(d).rstrip('0.') + "$", (0.95 * xmax, 0.5 * image.ymin() + 0.5 * image.ymax()), fontsize = 15, rgbcolor = (0, 0, 0)) return image #graph_of_sigma(2.0, 0.250, 0, 20).save('monic-0-20.png') #graphics_array([[graph_of_sigma(2.0, 0.6, 0, 100, True)], [graph_of_sigma(2.0, 0.3, 0, 100, True)], [graph_of_sigma(2.0, 0.1, 0, 100, True)]]).save('monic-0-100.png', figsize = [8, 12]) #graph_of_sigma(3.0, 0.5, 0, 50).save('monic-0-50.png') #graphics_array([[graph_of_sigma(1.0, 0.5, 0, 100, "lambda"), graph_of_sigma(1.0, 0.2, 0, 100, "lambda")], [graph_of_sigma(1.0, 0.1, 0, 100, "lambda"), graph_of_sigma(1.0, 0.05, 0, 100, "lambda")]]).save('monic-0-100_lambda.png', figsize = [12, 6]) #graphics_array([[graph_of_sigma(1.0, 0.75, 0, 100, "s"), graph_of_sigma(3.0, 0.75, 0, 100, "s")], [graph_of_sigma(6.0, 0.75, 0, 100, "s"), graph_of_sigma(9.0, 0.75, 0, 100, "s")]]).save('monic-0-100_s.png', figsize = [12, 6]) sig = sigma(2, 0.25) tbl = [] arg = [0.4 * i for i in range(0, 50)] for i in range(10): tbl += [[round(arg[i], 5), round(sig(arg[i]), 5), round(arg[10+i], 5), round(sig(arg[10+i]), 5), round(arg[20+i], 5), round(sig(arg[20+i]), 5), round(arg[30+i], 5), round(sig(arg[30+i]), 5), round(arg[40+i], 5), round(sig(arg[40+i]), 5)]] #latex(table(tbl, frame=True)) sig2 = sigma(3, 0.5) tbl2 = [] for i in range(10): tbl2 += [[round(i, 5), round(sig2(i), 5), round(10+i, 5), round(sig2(10+i), 5), round(20+i, 5), round(sig2(20+i), 5), round(30+i, 5), round(sig2(30+i), 5), round(40+i, 5), round(sig2(40+i), 5)]] #latex(table(tbl2, frame=True)) def print_table_of_parameters(f, a, b, errors): tbl = [] for epsilon in errors: c1, c2, theta1, theta2 = coefficients_of_approximation(1/4, f, a, b, epsilon) tbl += [[2, c1.n(25), c2.n(25), theta1.n(25), theta2, epsilon.n()]] print latex(table(tbl, header_row=True, frame=True)) #f(x) = 1 + x + x^2 / 2 + x^3 / 6 + x^4 / 24 + x^5 / 120 + x^6 / 720; print_table_of_parameters(f, -1, 1, [0.95, 0.6, 0.35, 0.1, 0.04, 0.01]) #f(x) = 4 * x / (4 + x^2); print_table_of_parameters(f, -1, 1, [0.95, 0.6, 0.35, 0.1, 0.04, 0.01]) #f(x) = sin(x) - x * cos(x + 1); print_table_of_parameters(f, -1, 1, [0.95, 0.6, 0.35, 0.1, 0.04, 0.01]) def graph_with_approximators(f, a, b, errors, colors, linestyles): image = plot(f, (x, a, b), color = 'blue', legend_label = '$f$') for i in range(len(errors)): p = polynomial_function(rational_approximating_Taylor_polynomial(f, a, b, errors[i])); pol = lambda t: p((t - a) / (b - a)) image += plot(pol, (x, a, b), color = colors[i], linestyle = linestyles[i], legend_label = '$\mathcal{N}\ $ ($\epsilon = ' + str(errors[i]).rstrip('0') + '$)') return image #f(x) = 1 + x + x^2 / 2 + x^3 / 6 + x^4 / 24 + x^5 / 120 + x^6 / 720; graph_with_approximators(f, -1, 1, [0.6, 0.35, 0.1], [(1, 0, 0), (2/3, 0, 1/3), (1/3, 0, 2/3)], [":", "-.", "--"]).save('monic-f1.png') #f(x) = 4 * x / (4 + x^2); graph_with_approximators(f, -1, 1, [0.6, 0.35, 0.1], [(1, 0, 0), (2/3, 0, 1/3), (1/3, 0, 2/3)], [":", "-.", "--"]).save('monic-f2.png') #f(x) = sin(x) - x * cos(x + 1); graph_with_approximators(f, -1, 1, [0.6, 0.35, 0.1], [(1, 0, 0), (2/3, 0, 1/3), (1/3, 0, 2/3)], [":", "-.", "--"]).save('monic-f3.png') @interact def _(a_b = range_slider(-10.0, 10.0, 1, (-1.0, 1.0), label = "[a, b]: "), lambda_ = slider(0.01, 1.0, step_size=0.001, default=0.2500, label="λ: "), x_range = range_slider(-10.0, 100.0, 1, (0, 20.0), label = "range: ")): d = a_b[1] - a_b[0] s = sigma(d, lambda_) h = lambda t: 1 - min(0.5, lambda_) / (log(t - d + 1) + 1) if t >= d else 1 - min(0.5, lambda_) xmin, xmax = x_range[0], x_range[1] show(plot(s, (x, xmin, xmax), color = 'blue') + plot(h, (x, xmin, xmax), color = 'red', linestyle = "--") + plot(1, (x, xmin, xmax), color = 'green', linestyle = "--")) #, aspect_ratio = 1 @interact def __(a_b = range_slider(-10.0, 10.0, 1, (-1.0, 1.0), label = "[a, b]: "), lambda_ = slider(0.01, 1.0, step_size=0.001, default=0.2500, label="λ: "), t = input_box(default = 0, label = "t: ")): d = a_b[1] - a_b[0] s = sigma(d, lambda_) sys.stdout.write('σ(') sys.stdout.write(str(t.n())) sys.stdout.write(') = ') sys.stdout.write(str(s(t).n())) @interact def ___(a_b = range_slider(-10.0, 10.0, 1, (-1.0, 1.0), label = "[a, b]: "), lambda_ = slider(0.01, 1.0, step_size=0.001, default=0.2500, label="λ: "), epsilon = slider(0.01, 1.0, step_size=0.01, default=1.0, label="ε: "), f = input_box(default = sin(x), label = "f: ")): p = rational_approximating_Taylor_polynomial(f, a_b[0], a_b[1], epsilon) if p == [0, 0]: show(plot(f, (x, a_b[0], a_b[1]), color = 'blue', legend_label = '$f$')) print "None of the first 100 Taylor polynomials approximates f with the given accuracy." else: polf = polynomial_function(p) pol1 = lambda t: polf((t - a_b[0]) / (a_b[1] - a_b[0])) show(plot(f, (x, a_b[0], a_b[1]), color = 'blue', legend_label = '$f$') + plot(pol1, (x, a_b[0], a_b[1]), color = 'red', linestyle = "--", legend_label = '$\mathcal{N}$')) c1, c2, theta1, theta2 = coefficients_of_approximation(lambda_, f, a_b[0], a_b[1], epsilon) if theta2 < 0: print "c₁ =", c1, ", c₂ =", c2, ", θ₁ =", theta1, ", θ₂ =", theta2 elif theta2 == 0: print "Error: Position of the polynomial is too large to calculate."
Interact: please open in CoCalc
Interact: please open in CoCalc
Interact: please open in CoCalc