[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
The fft
package comprises functions for the numerical (not symbolic)
computation of the fast Fourier transform.
Categories: Fourier transform · Numerical methods · Share packages · Package fft
[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
Translates complex values of the form r %e^(%i t)
to the form
a + b %i
, where r is the magnitude and t is the phase.
r and t are 1-dimensional arrays of the same size.
The array size need not be a power of 2.
The original values of the input arrays are
replaced by the real and imaginary parts, a
and b
, on return.
The outputs are calculated as
a = r cos(t) b = r sin(t)
polartorect
is the inverse function of recttopolar
.
load(fft)
loads this function. See also fft
.
Categories: Package fft · Complex variables
Translates complex values of the form a + b %i
to the form
r %e^(%i t)
, where a is the real part and b is the imaginary
part. a and b are 1-dimensional arrays of the same size.
The array size need not be a power of 2.
The original values of the input arrays are
replaced by the magnitude and angle, r
and t
, on return.
The outputs are calculated as
r = sqrt(a^2 + b^2) t = atan2(b, a)
The computed angle is in the range -%pi
to %pi
.
recttopolar
is the inverse function of polartorect
.
load(fft)
loads this function. See also fft
.
Categories: Package fft · Complex variables
Computes the inverse complex fast Fourier transform.
y is a list or array (named or unnamed) which contains the data to
transform. The number of elements must be a power of 2.
The elements must be literal numbers (integers, rationals, floats, or bigfloats)
or symbolic constants,
or expressions a + b*%i
where a
and b
are literal numbers
or symbolic constants.
inverse_fft
returns a new object of the same type as y,
which is not modified.
Results are always computed as floats
or expressions a + b*%i
where a
and b
are floats.
The inverse discrete Fourier transform is defined as follows.
Let x
be the output of the inverse transform.
Then for j
from 0 through n - 1
,
x[j] = sum(y[k] exp(-2 %i %pi j k / n), k, 0, n - 1)
load(fft)
loads this function.
See also fft
(forward transform), recttopolar
, and
polartorect
.
Examples:
Real data.
(%i1) load (fft) $ (%i2) fpprintprec : 4 $ (%i3) L : [1, 2, 3, 4, -1, -2, -3, -4] $ (%i4) L1 : inverse_fft (L); (%o4) [0.0, 14.49 %i - .8284, 0.0, 2.485 %i + 4.828, 0.0, 4.828 - 2.485 %i, 0.0, - 14.49 %i - .8284] (%i5) L2 : fft (L1); (%o5) [1.0, 2.0 - 2.168L-19 %i, 3.0 - 7.525L-20 %i, 4.0 - 4.256L-19 %i, - 1.0, 2.168L-19 %i - 2.0, 7.525L-20 %i - 3.0, 4.256L-19 %i - 4.0] (%i6) lmax (abs (L2 - L)); (%o6) 3.545L-16
Complex data.
(%i1) load (fft) $ (%i2) fpprintprec : 4 $ (%i3) L : [1, 1 + %i, 1 - %i, -1, -1, 1 - %i, 1 + %i, 1] $ (%i4) L1 : inverse_fft (L); (%o4) [4.0, 2.711L-19 %i + 4.0, 2.0 %i - 2.0, - 2.828 %i - 2.828, 0.0, 5.421L-20 %i + 4.0, - 2.0 %i - 2.0, 2.828 %i + 2.828] (%i5) L2 : fft (L1); (%o5) [4.066E-20 %i + 1.0, 1.0 %i + 1.0, 1.0 - 1.0 %i, 1.55L-19 %i - 1.0, - 4.066E-20 %i - 1.0, 1.0 - 1.0 %i, 1.0 %i + 1.0, 1.0 - 7.368L-20 %i] (%i6) lmax (abs (L2 - L)); (%o6) 6.841L-17
Categories: Package fft
Computes the complex fast Fourier transform.
x is a list or array (named or unnamed) which contains the data to
transform. The number of elements must be a power of 2.
The elements must be literal numbers (integers, rationals, floats, or bigfloats)
or symbolic constants,
or expressions a + b*%i
where a
and b
are literal numbers
or symbolic constants.
fft
returns a new object of the same type as x,
which is not modified.
Results are always computed as floats
or expressions a + b*%i
where a
and b
are floats.
The discrete Fourier transform is defined as follows.
Let y
be the output of the transform.
Then for k
from 0 through n - 1
,
y[k] = (1/n) sum(x[j] exp(+2 %i %pi j k / n), j, 0, n - 1)
When the data x are real,
real coefficients a
and b
can be computed such that
x[j] = sum(a[k]*cos(2*%pi*j*k/n)+b[k]*sin(2*%pi*j*k/n), k, 0, n/2)
with
a[0] = realpart (y[0]) b[0] = 0
and, for k from 1 through n/2 - 1,
a[k] = realpart (y[k] + y[n - k]) b[k] = imagpart (y[n - k] - y[k])
and
a[n/2] = realpart (y[n/2]) b[n/2] = 0
load(fft)
loads this function.
See also inverse_fft
(inverse transform), recttopolar
, and
polartorect
.
Examples:
Real data.
(%i1) load (fft) $ (%i2) fpprintprec : 4 $ (%i3) L : [1, 2, 3, 4, -1, -2, -3, -4] $ (%i4) L1 : fft (L); (%o4) [0.0, - 1.811 %i - .1036, 0.0, .6036 - .3107 %i, 0.0, .3107 %i + .6036, 0.0, 1.811 %i - .1036] (%i5) L2 : inverse_fft (L1); (%o5) [1.0, 2.168L-19 %i + 2.0, 7.525L-20 %i + 3.0, 4.256L-19 %i + 4.0, - 1.0, - 2.168L-19 %i - 2.0, - 7.525L-20 %i - 3.0, - 4.256L-19 %i - 4.0] (%i6) lmax (abs (L2 - L)); (%o6) 3.545L-16
Complex data.
(%i1) load (fft) $ (%i2) fpprintprec : 4 $ (%i3) L : [1, 1 + %i, 1 - %i, -1, -1, 1 - %i, 1 + %i, 1] $ (%i4) L1 : fft (L); (%o4) [0.5, .3536 %i + .3536, - 0.25 %i - 0.25, 0.5 - 6.776L-21 %i, 0.0, - .3536 %i - .3536, 0.25 %i - 0.25, 0.5 - 3.388L-20 %i] (%i5) L2 : inverse_fft (L1); (%o5) [1.0 - 4.066E-20 %i, 1.0 %i + 1.0, 1.0 - 1.0 %i, - 1.008L-19 %i - 1.0, 4.066E-20 %i - 1.0, 1.0 - 1.0 %i, 1.0 %i + 1.0, 1.947L-20 %i + 1.0] (%i6) lmax (abs (L2 - L)); (%o6) 6.83L-17
Computation of sine and cosine coefficients.
(%i1) load (fft) $ (%i2) fpprintprec : 4 $ (%i3) L : [1, 2, 3, 4, 5, 6, 7, 8] $ (%i4) n : length (L) $ (%i5) x : make_array (any, n) $ (%i6) fillarray (x, L) $ (%i7) y : fft (x) $ (%i8) a : make_array (any, n/2 + 1) $ (%i9) b : make_array (any, n/2 + 1) $ (%i10) a[0] : realpart (y[0]) $ (%i11) b[0] : 0 $ (%i12) for k : 1 thru n/2 - 1 do (a[k] : realpart (y[k] + y[n - k]), b[k] : imagpart (y[n - k] - y[k])); (%o12) done (%i13) a[n/2] : y[n/2] $ (%i14) b[n/2] : 0 $ (%i15) listarray (a); (%o15) [4.5, - 1.0, - 1.0, - 1.0, - 0.5] (%i16) listarray (b); (%o16) [0, - 2.414, - 1.0, - .4142, 0] (%i17) f(j) := sum (a[k]*cos(2*%pi*j*k/n) + b[k]*sin(2*%pi*j*k/n), k, 0, n/2) $ (%i18) makelist (float (f (j)), j, 0, n - 1); (%o18) [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
Categories: Package fft
[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
Returns a rearranged representation of expr as in Horner's rule, using
x as the main variable if it is specified. x
may be omitted in
which case the main variable of the canonical rational expression form of
expr is used.
horner
sometimes improves stability if expr
is
to be numerically evaluated. It is also useful if Maxima is used to
generate programs to be run in Fortran. See also stringout
.
(%i1) expr: 1e-155*x^2 - 5.5*x + 5.2e155; 2 (%o1) 1.e-155 x - 5.5 x + 5.2e+155 (%i2) expr2: horner (%, x), keepfloat: true; (%o2) 1.0 ((1.e-155 x - 5.5) x + 5.2e+155) (%i3) ev (expr, x=1e155); Maxima encountered a Lisp error: arithmetic error FLOATING-POINT-OVERFLOW signalled Automatically continuing. To enable the Lisp debugger set *debugger-hook* to nil. (%i4) ev (expr2, x=1e155); (%o4) 7.00000000000001e+154
Categories: Numerical methods
Finds a root of the expression expr or the function f over the
closed interval [a, b]. The expression expr may be an
equation, in which case find_root
seeks a root of
lhs(expr) - rhs(expr)
.
Given that Maxima can evaluate expr or f over
[a, b] and that expr or f is continuous,
find_root
is guaranteed to find the root,
or one of the roots if there is more than one.
find_root
initially applies binary search.
If the function in question appears to be smooth enough,
find_root
applies linear interpolation instead.
bf_find_root
is a bigfloat version of find_root
. The
function is computed using bigfloat arithmetic and a bigfloat result
is returned. Otherwise, bf_find_root
is identical to
find_root
, and the following description is equally applicable
to bf_find_root
.
The accuracy of find_root
is governed by abserr
and
relerr
, which are optional keyword arguments to
find_root
. These keyword arguments take the form
key=val
. The keyword arguments are
Desired absolute error of function value at root. Default is
find_root_abs
.
Desired relative error of root. Default is find_root_rel
.
find_root
stops when the function in question evaluates to
something less than or equal to abserr
, or if successive
approximants x_0, x_1 differ by no more than relerr
* max(abs(x_0), abs(x_1))
. The default values of
find_root_abs
and find_root_rel
are both zero.
find_root
expects the function in question to have a different sign at
the endpoints of the search interval.
When the function evaluates to a number at both endpoints
and these numbers have the same sign,
the behavior of find_root
is governed by find_root_error
.
When find_root_error
is true
,
find_root
prints an error message.
Otherwise find_root
returns the value of find_root_error
.
The default value of find_root_error
is true
.
If f evaluates to something other than a number at any step in the search
algorithm, find_root
returns a partially-evaluated find_root
expression.
The order of a and b is ignored; the region in which a root is sought is [min(a, b), max(a, b)].
Examples:
(%i1) f(x) := sin(x) - x/2; x (%o1) f(x) := sin(x) - - 2 (%i2) find_root (sin(x) - x/2, x, 0.1, %pi); (%o2) 1.895494267033981 (%i3) find_root (sin(x) = x/2, x, 0.1, %pi); (%o3) 1.895494267033981 (%i4) find_root (f(x), x, 0.1, %pi); (%o4) 1.895494267033981 (%i5) find_root (f, 0.1, %pi); (%o5) 1.895494267033981 (%i6) find_root (exp(x) = y, x, 0, 100); x (%o6) find_root(%e = y, x, 0.0, 100.0) (%i7) find_root (exp(x) = y, x, 0, 100), y = 10; (%o7) 2.302585092994046 (%i8) log (10.0); (%o8) 2.302585092994046 (%i9) fpprec:32; (%o9) 32 (%i10) bf_find_root (exp(x) = y, x, 0, 100), y = 10; (%o10) 2.3025850929940456840179914546844b0 (%i11) log(10b0); (%o11) 2.3025850929940456840179914546844b0
Categories: Algebraic equations · Numerical methods
Returns an approximate solution of expr = 0
by Newton's method,
considering expr to be a function of one variable, x.
The search begins with x = x_0
and proceeds until abs(expr) < eps
(with expr evaluated at the current value of x).
newton
allows undefined variables to appear in expr,
so long as the termination test abs(expr) < eps
evaluates
to true
or false
.
Thus it is not necessary that expr evaluate to a number.
load(newton1)
loads this function.
See also realroots
, allroots
, find_root
, and
mnewton
.
Examples:
(%i1) load (newton1); (%o1) /maxima/share/numeric/newton1.mac (%i2) newton (cos (u), u, 1, 1/100); (%o2) 1.570675277161251 (%i3) ev (cos (u), u = %); (%o3) 1.2104963335033529e-4 (%i4) assume (a > 0); (%o4) [a > 0] (%i5) newton (x^2 - a^2, x, a/2, a^2/100); (%o5) 1.00030487804878 a (%i6) ev (x^2 - a^2, x = %); 2 (%o6) 6.098490481853958e-4 a
Categories: Algebraic equations · Numerical methods
[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
The Ordinary Differential Equations (ODE) solved by the functions in this section should have the form,
dy -- = F(x,y) dx
which is a first-order ODE. Higher order differential equations of order n must be written as a system of n first-order equations of that kind. For instance, a second-order ODE should be written as a system of two equations
dx dy -- = G(x,y,t) -- = F(x,y,t) dt dt
The first argument in the functions will be a list with the expressions on the right-side of the ODE's. The variables whose derivatives are represented by those expressions should be given in a second list. In the case above those variables are x and y. The independent variable, t in the examples above, might be given in a separated option. If the expressions given do not depend on that independent variable, the system is called autonomous.
Categories: Differential equations · Numerical methods · Plotting
[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
[
u,v]
, ...options...)
[
dxdt,dydt]
, ...options...)
[
dudt,dvdt]
, [
u,v]
, ...options...)
The function plotdf
creates a two-dimensional plot of the direction
field (also called slope field) for a first-order Ordinary Differential
Equation (ODE) or a system of two autonomous first-order ODE's.
Plotdf requires Xmaxima, even if its run from a Maxima session in a console, since the plot will be created by the Tk scripts in Xmaxima. If Xmaxima is not installed plotdf will not work.
dydx, dxdt and dydt are expressions that depend on
x and y. dvdu, dudt and dvdt are
expressions that depend on u and v. In addition to those two
variables, the expressions can also depend on a set of parameters, with
numerical values given with the parameters
option (the option
syntax is given below), or with a range of allowed values specified by a
sliders option.
Several other options can be given within the command, or selected in
the menu. Integral curves can be obtained by clicking on the plot, or
with the option trajectory_at
. The direction of the integration
can be controlled with the direction
option, which can have
values of forward, backward or both. The number of
integration steps is given by nsteps
and the time interval
between them is set up with the tstep
option. The Adams Moulton
method is used for the integration; it is also possible to switch to an
adaptive Runge-Kutta 4th order method.
Plot window menu:
The menu in the plot window has the following options: Zoom, will change the behavior of the mouse so that it will allow you to zoom in on a region of the plot by clicking with the left button. Each click near a point magnifies the plot, keeping the center at the point where you clicked. Holding the Shift key while clicking, zooms out to the previous magnification. To resume computing trajectories when you click on a point, select Integrate from the menu.
The option Config in the menu can be used to change the ODE(s) in use and various other settings. After configuration changes are made, the menu option Replot should be selected, to activate the new settings. If a pair of coordinates are entered in the field Trajectory at in the Config dialog menu, and the enter key is pressed, a new integral curve will be shown, in addition to the ones already shown. When Replot is selected, only the last integral curve entered will be shown.
Holding the right mouse button down while the cursor is moved, can be used to drag the plot sideways or up and down. Additional parameters such as the number of steps, the initial value of t and the x and y centers and radii, may be set in the Config menu.
A copy of the plot can be saved as a postscript file, using the menu option Save.
Plot options:
The plotdf
command may include several options, each one being
a list of two or more elements. The first element in each option is the name
of the option, and the remainder is the value or values assigned to the
option.
The options which are recognized by plotdf
are the following:
plotdf
, the x
variable will be directly proportional to t.
The default value is 0.1.
tstep
that will be used for the independent variable, to compute an integral
curve.
The default value is 100.
forward
, to make the independent variable increase
nsteps
times, with increments tstep
, backward
, to
make the independent variable decrease, or both
that will lead to
an integral curve that extends nsteps
forward, and nsteps
backward. The keywords right
and left
can be used as
synonyms for forward
and backward
.
The default value is both
.
versus_t
is given any value
different from 0, the second plot window will be displayed. The second
plot window includes another menu, similar to the menu of the main plot
window.
The default value is 0.
name=value
.
name=min:max
Examples:
(%i1) plotdf(exp(-x)+y,[trajectory_at,2,-0.1])$
(%i1) plotdf(x-y^2,[xfun,"sqrt(x);-sqrt(x)"], [trajectory_at,-1,3], [direction,forward], [y,-5,5], [x,-4,16])$
The graph also shows the function y = sqrt(x).
(%i1) plotdf([v,-k*z/m], [z,v], [parameters,"m=2,k=2"], [sliders,"m=1:5"], [trajectory_at,6,0])$
(%i1) plotdf([y,-(k*x + c*y + b*x^3)/m], [parameters,"k=-1,m=1.0,c=0,b=1"], [sliders,"k=-2:2,m=-1:1"],[tstep,0.1])$
(%i1) plotdf([w,-g*sin(a)/l - b*w/m/l], [a,w], [parameters,"g=9.8,l=0.5,m=0.3,b=0.05"], [trajectory_at,1.05,-9],[tstep,0.01], [a,-10,2], [w,-14,14], [direction,forward], [nsteps,300], [sliders,"m=0.1:1"], [versus_t,1])$
Categories: Differential equations · Plotting · Numerical methods
Plots equipotential curves for exp, which should be an expression depending on two variables. The curves are obtained by integrating the differential equation that define the ortogonal trajetories to the solutions of the autonomous system obtained from the gradient of the expression given. The plot can also show the integral curves for that gradient system (option fieldlines).
This program also requires Xmaxima, even if its run from a Maxima session in a console, since the plot will be created by the Tk scripts in Xmaxima. By default, the plot region will be empty until the user clicks in a point (or gives its coordinate with in the set-up menu or via the trajectory_at option).
Most options accepted by plotdf can also be used for ploteq and the plot interface is the same that was described in plotdf.
Example:
(%i1) V: 900/((x+1)^2+y^2)^(1/2)-900/((x-1)^2+y^2)^(1/2)$ (%i2) ploteq(V,[x,-2,2],[y,-2,2],[fieldlines,"blue"])$
clicking on a point will plot the equipotential curve that passes by that point (in red) and the ortogonal trajectory (in blue).
Categories: Differential equations · Plotting · Numerical methods
The first form solves numerically one first-order ordinary differential equation, and the second form solves a system of m of those equations, using the 4th order Runge-Kutta method. var represents the dependent variable. ODE must be an expression that depends only on the independent and dependent variables and defines the derivative of the dependent variable with respect to the independent variable.
The independent variable is specified with domain
, which must be a
list of four elements as, for instance:
[t, 0, 10, 0.1]
the first element of the list identifies the independent variable, the second and third elements are the initial and final values for that variable, and the last element sets the increments that should be used within that interval.
If m equations are going to be solved, there should be m
dependent variables v1, v2, ..., vm. The initial values
for those variables will be init1, init2, ..., initm.
There will still be just one independent variable defined by domain
,
as in the previous case. ODE1, ..., ODEm are the expressions
that define the derivatives of each dependent variable in
terms of the independent variable. The only variables that may appear in
those expressions are the independent variable and any of the dependent
variables. It is important to give the derivatives ODE1, ...,
ODEm in the list in exactly the same order used for the dependent
variables; for instance, the third element in the list will be interpreted
as the derivative of the third dependent variable.
The program will try to integrate the equations from the initial value of the independent variable until its last value, using constant increments. If at some step one of the dependent variables takes an absolute value too large, the integration will be interrupted at that point. The result will be a list with as many elements as the number of iterations made. Each element in the results list is itself another list with m+1 elements: the value of the independent variable, followed by the values of the dependent variables corresponding to that point.
To solve numerically the differential equation
dx/dt = t - x^2
With initial value x(t=0) = 1, in the interval of t from 0 to 8 and with increments of 0.1 for t, use:
(%i1) results: rk(t-x^2,x,1,[t,0,8,0.1])$ (%i2) plot2d ([discrete, results])$
the results will be saved in the list results
and the plot will show the solution obtained, with t on the horizontal axis and x on the vertical axis.
To solve numerically the system:
dx/dt = 4-x^2-4*y^2 dy/dt = y^2-x^2+1
for t between 0 and 4, and with values of -1.25 and 0.75 for x and y at t=0:
(%i1) sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$ (%i2) plot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$
The plot will show the solution for variable y as a function of t.
Categories: Differential equations · Numerical methods
[ << ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
This document was generated by Jaime Villate on April, 11 2013 using texi2html 1.76.