---
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
-~----------~----~----~----~------~----~------~--~---
