--- examples/all | 1 + examples/gibbs_phenomenon.py | 131 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 132 insertions(+), 0 deletions(-) create mode 100755 examples/gibbs_phenomenon.py
diff --git a/examples/all b/examples/all index f2ad056..43ac2dc 100755 --- a/examples/all +++ b/examples/all @@ -23,6 +23,7 @@ ./differential_equations.py #./pidigits.py ./differentiation.py +./gibbs_phenomenon.py #./plotting_nice_plot.py ./series.py ./expansion.py diff --git a/examples/gibbs_phenomenon.py b/examples/gibbs_phenomenon.py new file mode 100755 index 0000000..64f8c02 --- /dev/null +++ b/examples/gibbs_phenomenon.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python + +""" +This example illustrates the Gibbs phenomenon. + +It also calculates the Wilbraham-Gibbs constant by two approaches: + +1) calculating the fourier series of the step function and determining the +first maximum. +2) evaluating the integral for si(pi). +""" + +import iam_sympy_example + +from sympy import var, sqrt, integrate, conjugate, seterr, abs, pprint, I, pi,\ + sin, cos, sign, Plot, msolve, lambdify, Integral, S + +#seterr(True) + +x = var("x", real=True) + +def l2_norm(f, lim): + """ + Calculates L2 norm of the function "f", over the domain lim=(x, a, b). + + x ...... the independent variable in f over which to integrate + a, b ... the limits of the interval + + Example: + + >>> l2_norm(1, (x, -1, 1)) + 2**(1/2) + >>> l2_norm(x, (x, -1, 1)) + 1/3*6**(1/2) + """ + return sqrt(integrate(abs(f)**2, lim)) + +def l2_inner_product(a, b, lim): + """ + Calculates the L2 inner product (a, b) over the domain lim. + """ + return integrate(conjugate(a)*b, lim) + +def l2_projection(f, basis, lim): + """ + L2 projects the function f on the basis over the domain lim. + """ + r = 0 + for b in basis: + r += l2_inner_product(f, b, lim) * b + return r + +def l2_gram_schmidt(list, lim): + """ + Orthonormalizes the "list" of functions using the Gram-Schmidt process. + + Example: + >>> l2_gram_schmidt([1, x, x**2], (x, -1, 1)] + [1/2*2**(1/2), x*6**(1/2)/2, -3*10**(1/2)*(1/3 - x**2)/4] + """ + r = [] + for a in list: + if r == []: + v = a + else: + v = a - l2_projection(a, r, lim) + v_norm = l2_norm(v, lim) + if v_norm == 0: + raise Exception("The sequence is not linearly independent.") + r.append(v/v_norm) + return r + +#l = l2_gram_schmidt([1, cos(x), sin(x), cos(2*x), sin(2*x)], (x, -pi, pi)) +#l = l2_gram_schmidt([1, cos(x), sin(x)], (x, -pi, pi)) +# the code below is equivalen to l2_gram_schmidt(), but faster: +l = [1/sqrt(2)] +for i in range(1, 100): + l.append(cos(i*x)) + l.append(sin(i*x)) +l = [f/sqrt(pi) for f in l] + +def integ(f): + return integrate(f, (x, -pi, 0)) + integrate(-f, (x, 0, pi)) + +def series(l): + """ + Normalizes the series. + """ + r = 0 + for b in l: + r += integ(b)*b + return r + +def msolve(f, x): + """ + Finds the first root of f(x) to the left of 0. + + The x0 and dx below are taylored to get the correct result for our + particular function --- the general solver often overshoots the first + solution. + """ + f = lambdify(x, f) + x0 = -0.001 + dx = 0.001 + while f(x0-dx) * f(x0) > 0: + x0 = x0-dx + x_max = x0-dx + x_min = x0 + assert f(x_max) > 0 + assert f(x_min) < 0 + for n in range(100): + x0 = (x_max+x_min)/2 + if f(x0) > 0: + x_max = x0 + else: + x_min = x0 + return x0 + +f = series(l) +print "Fourier series of the step function" +pprint(f) +#Plot(f.diff(x), [x, -5, 5, 3000]) +x0 = msolve(f.diff(x), x) + +print "x-value of the maximum:", x0 +max = f.subs(x, x0).evalf() +print "y-value of the maximum:", max +g = max*pi/2 +print "Wilbraham-Gibbs constant :", g.evalf() +print "Wilbraham-Gibbs constant (exact):", \ + Integral(sin(x)/x, (x, 0, pi)).evalf() -- 1.5.6.5 --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy-patches" group. To post to this group, send email to sympy-patches@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sympy-patches?hl=en -~----------~----~----~----~------~----~------~--~---