Feeds:
Posts

## SAGE tip: Fourier Series Approximation with Interactive mode

After my last post, about Fouries Series, I decided to give a try to the interactive mode of SAGE.

🙂

What did I do?

First I check out the examples worked at the Wiki page. In the section of graphics I got inspired from the post Interactive 2D Plotting by Timothy Clemans and Newton’s Method by William Stein post in the Calculus section. So… there it is!

The code is too long to be posted here… but here it goes anyway

var('x,n,i')

def error_msg(msg): ... print '<p style="font-family: Arial, sans-serif; color: #000;"><span     style="color: red; font-weight: bold;">Error</span>: %s</p>' % msg

...

def a(f, n, init, final): ... x = f.variables()[0] ... L = (final - init)/2 ... coeff = integral(f(x)*cos(n*pi*x/L), (x,init,final))/L ... return coeff

 ... def b(f, n, init, final): ... x = f.variables()[0] ... L = (final - init)/2 ... coeff = integral(f(x)*sin(n*pi*x/L), (x,init,final))/L ... return coeff ... def FS(f, n, init, final): ... x = f.variables()[0] ... L = (final - init)/2 ... return a(f, 0, init, final)/2 + sum(a(f, i, init, final)*cos(i*x*pi/L)+b(f, i, init, final)*sin(i*x*pi/L), i, 1, n) ... 

@interact def interactive_2d_plotter(clr=Color('green'), expression=input_box('x^2', 'Expression', str), x_range=range_slider(-6,6,1,(-4,4), label='X Range'), order=slider([0..100],None,None,3), square=checkbox(True, 'Square'), axes=checkbox(False, 'Show Axes')): ... if expression: ...   try: ...      expression = SR(expression) # turn string into a Sage expression ...   except TypeError: ...      print error_msg('This is not an expression.') ...      return ...   try: ...      xmin, xmax = x_range ...      n = order ...      if square or not axes: ...         print "var('%s')\nplot(%s).show(%s%s%s)" % (expression.variables()[0], repr(expression), 'aspect_ratio=1' if square else '', ', ' if square and not axes else '', 'axes=False' if not axes else '') ...         if square: ...            P = plot(FS(expression, n, xmin, xmax), xmin, xmax, color=clr) + plot(expression, xmin, xmax) ...            P.show(aspect_ratio=1, axes=axes) ...         else: ...            P = plot(FS(expression, n, xmin, xmax), xmin, xmax, color=clr) + plot(expression, xmin, xmax) ...            P.show(axes=axes) ...      else: ...         print "var('%s')\nplot(%s)" % (expression.variables()[0], repr(expression)) ...         P = plot(FS(expression, n, xmin, xmax), xmin, xmax, color=clr) + plot(expression, xmin, xmax) ...         P.show(axes=axes) ...   except ValueError: ...      print error_msg('This expression has more than one variable.') ...      return ...   except TypeError: ...      print error_msg("This expression contains an unknown function.") ...      return

because it’s almost in-understandable 😛

Some pictures!!!

The code is set to give the Fourier approximation of the $f(x)=x^2$ with aspect_ratio = (1,1).

Initial Plot. Approx. of x^2

One can set a non one-to-one aspect_ratio, by un-checking the box,

Initial Plot with-no-square

And we can view the frame as well by checking the other box,

Initial plot no-square and frame

The first box allows us to change the colour, by using a palette,

Initial plot with color-palette

for example, set the colour to red,

Initial plot in red

We can also change the function, to an exponential, $f(x)=e^x$,

Fourier approximation of the exponential

and change the domain of the function,

Changing the interval

or the order of the approximation,

Exponential changing the order of the approximation

I published my sage file on my page (at the end of the page in the attached files)… for the sake of clarity (and completeness).

Enjoy,

Dox

## SAGE tip: Fourier Series Approximation

Inspired by a post in sage-devel (or support) group of SAGE, I came along with this few lines which allows me to plot a Fourier Series Approximation of the line, to a given order,

sage: reset() sage: var('x,i,n') (x,i,n) sage: def b(n): ... return 2.*(-1)^(n+1)/n ... sage: def f(x,n): ... return sum(b(i)*sin(i*x),i,1,n) ... sage: p = Graphics() sage: for n in [1,2,3,4,5,6]: ... p = plot(f(x,n), (x,-pi,pi), color=hue((n+1)/7.0)) + plot(x, -pi, pi, color='black') + text("Fourier Approximation of order %d" %n, (-1,3), fontsize=14, color='blue') ... p ...

and get as result, the series of plots that follows,

1st Fourier Approximation to a Line

2nd Fourier Approximation to a Line

3rd Fourier Approximation to a Line

4th Fourier Approximation to a Line

5th Fourier Approximation to a Line

6th Fourier Approximation to a Line

However, I was not happy ’cause I introduce the Fourier coefficients… So, I did a second try.

sage: reset() sage: var('x,n,i') (x,n,i) sage: f(x) = x^2

 sage: def a(n): ... coeff = integral(f(x)*cos(n*x), (x,-pi,pi))/pi ... return coeff ... sage: def b(n): ... coeff = integral(f(x)*sin(n*x), (x,-pi,pi))/pi ... return coeff ... sage: def FS(n): ... return a(0)/2 + sum(a(i)*cos(i*x)+b(i)*sin(i*x), i, 1, n) ... 

Where I’ve defined the Fourier coefficients and the Fourier Series of a given function $f(x)$, introduced by the user.

With the line

sage: for n in [1,2,3,4,5,6]: ... p = plot(FS(n), (x,-1.1*pi,1.1*pi), color=hue((n+1)/7.0)) + plot(f(x), (x, -1.1*pi, 1.1*pi), color='black') + text("Fourier Approximation of order %d" %n, (0,-2), fontsize=14, color='blue') ... p ...

One get the Fourier Series Approximation for $f(x)=x^2$,

1st Approximation for the parabola

2nd Approximation for the parabola

3rd Approximation for the parabola

4th Approximation for the parabola

5th Approximation for the parabola

6th Approximation for the parabola

#### NOTE…

• The use of the last program is restricted to the interval $[-\pi,\pi]$.
• In order to find the Fourier coefficients the integrals might be doable… so no every function $f(x)$ can be shown as a Fourier series approximation.
• I tried to define the above using numerical_integrategral command, but didn’t work. Does anyone knows why?
• I also tried to use range command instead of a list for the loop… Didn’t work!!! Any clues?
• Ok, that was it. Enjoy!!!

DOX.

This morning I read a post at sage-devel group titled “Integral Functions”, and was really interesting.

For example, say you want to define a function as the primitive integral of another function,
$F(x) = \int_0^x\;f(x')\;dx'$,

sage: var('x') sage: f = function('f',x) sage: F = integral(f, x, 0, x) sage: F.derivative(x)

gives an error… Why? derivative expect all parameters of $F$ to be symbols, and zero is an integer.

## Bypass

A sort of solution to this problem is giving a whole set of parameters to the limits, say, $F(x, a, b) = \int_a^b\;dx\; f(x)$

sage: var('x,a,b') (x, a, b) sage: f = function('f',x) sage: F = integral(f, x, a, b) sage: F.derivative(x) 0 sage: F.derivative(a) -f(a) sage: F.derivative(b) f(b)

Thus, it’s well defined for any case.

If you don’t give integration limits,

sage: var('x') sage: f = function('f',x) sage: F = integral(f, x) sage: F.derivative(x) f(x)

it work as expected.

The discussion page is HERE

Enjoy!

DOX

## SAGE tip: rewriting expressions

Today I was calculating some stuff with the help of SAGE, and realize that the expressions got a lot (really, a lot) simpler if they where written in terms of  hyperbolic functions instead of exponentials.

$e^x = \cosh(x)+\sinh(x)$

$e^{-x} = \cosh(x)-\sinh(x)$

So I enter the sage-devel channel of IRC (freenode), but there was a lazy day around… Sunday. Therefore I decided to write to the sage-dev group on Google groups.

Francois Maltey answer my question on how to do the transformation… He has written a package that does it!

• Go to page http://wiki.sagemath.org/symbolics/rewrite
look at the second line the “attachement” (in smaller characters)
and get the most recent file.
• Then put this file in your Sage directory
In Sage, type : load “/the/good/file/in/the/good/directory’
Then call the rewrite function.
• Thank you Francois

#### Example

Suppose the path to the file is “/home/me/rewrite-xxx.sage”

sage: load "/home/me/rewrite-xxx.sage" sage: A = exp(x) + exp(-x) sage: rewrite(A, 'exp2sinhconh') 2*cosh(x)

If you go to http://wiki.sagemath.org/symbolics/rewrite, will find all possible commands which perform such transformations.

Enjoy!

Dox