[PATCH] gibbs_phenomenon example added (#1114)
Ondrej Certik
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
Andy Ray Terrel
different algorithms used for examples. It seems that we could really
just have such a large list of examples it is hard to keep up with
them. Really the idea I would have is to make a corresponding wiki
that explained the examples in more detail for someone trying to learn
things.
For example if I didn't know what the gram-schmidt process did, I
would be utterly confused as why your trick worked and what it was
doing.
Because you just give the basis later function l2_norm,
l2_inner_product, l2_projection, and l2_gram_schmidt are never used
and should be, IMHO, removed.
Ondrej Certik
>
> One thing that would be nice is to give a bit more information on
> different algorithms used for examples. It seems that we could really
> just have such a large list of examples it is hard to keep up with
> them. Really the idea I would have is to make a corresponding wiki
> that explained the examples in more detail for someone trying to learn
> things.
>
> For example if I didn't know what the gram-schmidt process did, I
> would be utterly confused as why your trick worked and what it was
> doing.
Right. I created an issue for that:
http://code.google.com/p/sympy/issues/detail?id=1196
>
> Because you just give the basis later function l2_norm,
> l2_inner_product, l2_projection, and l2_gram_schmidt are never used
> and should be, IMHO, removed.
Well, it's an example, so I wanted to leave it there so that people
know how it works. But it's true that it is never used.
Because if "f" is complex, you either need to do: f * f.conjugate()
or abs(f)**2.
Well, shouldn't rather one be using fonts that distinguish l, 1, i,
etc.? I think it is absolutely essential to being able to distinguish
between those.
Ondrej
Ondrej Certik
Andy --- are you ok with pushing this patch in and possibly improve
further, or should this be first improved and only then pushed in?
Ondrej
Vinzent Steinberg
background information. Some links would be enough.
Vinzent
On 14 Nov., 15:17, "Ondrej Certik" <ond...@certik.cz> wrote:
> On Fri, Nov 14, 2008 at 1:48 PM, Ondrej Certik <ond...@certik.cz> wrote: